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*/