Actual source code: ex3.c
1: static char help[] = "Tests LETKF-specific functionality: localization radius and type set/get.\n\n";
2: #include <petscda.h>
4: int main(int argc, char **argv)
5: {
6: PetscDA da;
7: PetscReal radius = 3.14, radius_check;
8: PetscDALETKFLocalizationType type_check;
10: PetscFunctionBeginUser;
11: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
13: /* Create the LETKF DA object */
14: PetscCall(PetscDACreate(PETSC_COMM_WORLD, &da));
15: PetscCall(PetscDASetSizes(da, 10, 10));
16: PetscCall(PetscDAEnsembleSetSize(da, 5));
17: PetscCall(PetscDASetFromOptions(da));
18: PetscCall(PetscDASetUp(da));
20: /* Test localization radius set/get round-trip. Exact == compare is intentional:
21: the setter stores the value verbatim and the getter returns it with no arithmetic. */
22: PetscCall(PetscDALETKFSetLocalizationRadius(da, radius));
23: PetscCall(PetscDALETKFGetLocalizationRadius(da, &radius_check));
24: PetscCheck(radius_check == radius, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "SetLocalizationRadius/GetLocalizationRadius round-trip failed: set %g, got %g", (double)radius, (double)radius_check);
25: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Localization radius set/get: %g\n", (double)radius_check));
27: /* Test localization type set/get round-trip across every enum value, and verify that
28: changing the type leaves the previously-set radius untouched. Loop in PetscInt to keep
29: compilers that treat the example as C++ (e.g. Apple clang) from rejecting enum++. */
30: for (PetscInt ti = PETSCDA_LETKF_LOC_NONE; ti < PETSCDA_LETKF_LOC_NUM_TYPES; ++ti) {
31: PetscDALETKFLocalizationType t = (PetscDALETKFLocalizationType)ti;
33: PetscCall(PetscDALETKFSetLocalizationType(da, t));
34: PetscCall(PetscDALETKFGetLocalizationType(da, &type_check));
35: PetscCheck(type_check == t, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "SetLocalizationType/GetLocalizationType round-trip failed at type %d: got %d", (int)t, (int)type_check);
36: PetscCall(PetscDALETKFGetLocalizationRadius(da, &radius_check));
37: PetscCheck(radius_check == radius, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Setting localization type %d clobbered radius: expected %g, got %g", (int)t, (double)radius, (double)radius_check);
38: }
39: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Localization type set/get round-trip: ok\n"));
41: PetscCall(PetscDADestroy(&da));
42: PetscCall(PetscFinalize());
43: return 0;
44: }
46: /*TEST
48: test:
49: suffix: 1
50: requires: !complex
52: TEST*/