Actual source code: petscda.c

      1: #include <petsc/private/daimpl.h>
      2: #include <petscblaslapack.h>

      4: PetscClassId      PETSCDA_CLASSID          = 0;
      5: PetscLogEvent     PetscDA_Analysis         = 0;
      6: PetscBool         PetscDARegisterAllCalled = PETSC_FALSE;
      7: PetscFunctionList PetscDAList              = NULL;

      9: static PetscBool PetscDAPackageInitialized = PETSC_FALSE;

     11: /*@C
     12:   PetscDAInitializePackage - This function initializes everything in the `PetscDA`
     13:   package. called on the first call to `PetscDACreate()` when using static or shared
     14:   libraries.

     16:   Logically Collective

     18:   Level: developer

     20: .seealso: `PetscDAFinalizePackage()`, `PetscInitialize()`
     21: @*/
     22: PetscErrorCode PetscDAInitializePackage(void)
     23: {
     24:   PetscFunctionBegin;
     25:   if (PetscDAPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);

     27:   PetscDAPackageInitialized = PETSC_TRUE;
     28:   PetscCall(PetscClassIdRegister("Data Assimilation", &PETSCDA_CLASSID));
     29:   PetscCall(PetscDARegisterAll());
     30:   PetscCall(PetscRegisterFinalize(PetscDAFinalizePackage));
     31:   PetscCall(PetscLogEventRegister("PetscDAAnalysis", PETSCDA_CLASSID, &PetscDA_Analysis));
     32:   PetscFunctionReturn(PETSC_SUCCESS);
     33: }

     35: /*@C
     36:   PetscDAFinalizePackage - This function finalizes everything in the `PetscDA` package. It
     37:   is called from `PetscFinalize()`.

     39:   Logically Collective

     41:   Level: developer

     43: .seealso: `PetscDAInitializePackage()`, `PetscInitialize()`
     44: @*/
     45: PetscErrorCode PetscDAFinalizePackage(void)
     46: {
     47:   PetscFunctionBegin;
     48:   PetscCall(PetscFunctionListDestroy(&PetscDAList));
     49:   PetscDARegisterAllCalled  = PETSC_FALSE;
     50:   PetscDAPackageInitialized = PETSC_FALSE;
     51:   PetscDA_Analysis          = 0;
     52:   PetscFunctionReturn(PETSC_SUCCESS);
     53: }

     55: /*@C
     56:   PetscDARegister - Registers a constructor for a `PetscDA` implementation with the
     57:   dispatcher.

     59:   Not Collective

     61:   Input Parameters:
     62: + sname    - name associated with the implementation
     63: - function - routine that creates the implementation and installs method table

     65:   Level: developer

     67: .seealso: [](ch_da), `PetscDARegisterAll()`, `PetscDASetType()`
     68: @*/
     69: PetscErrorCode PetscDARegister(const char sname[], PetscErrorCode (*function)(PetscDA))
     70: {
     71:   PetscFunctionBegin;
     72:   PetscCall(PetscDAInitializePackage());
     73:   PetscCall(PetscFunctionListAdd(&PetscDAList, sname, function));
     74:   PetscFunctionReturn(PETSC_SUCCESS);
     75: }

     77: PETSC_INTERN PetscErrorCode PetscDACreate_ETKF(PetscDA);
     78: PETSC_INTERN PetscErrorCode PetscDACreate_LETKF(PetscDA);

     80: /*@
     81:   PetscDARegisterAll - Registers all data assimilation backends that were compiled in.

     83:   Not Collective

     85:   Level: developer

     87: .seealso: [](ch_da), `PetscDARegister()`
     88: @*/
     89: PetscErrorCode PetscDARegisterAll(void)
     90: {
     91:   PetscFunctionBegin;
     92:   if (PetscDARegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
     93:   PetscDARegisterAllCalled = PETSC_TRUE;
     94:   PetscCall(PetscDARegister(PETSCDAETKF, PetscDACreate_ETKF));
     95:   PetscCall(PetscDARegister(PETSCDALETKF, PetscDACreate_LETKF));
     96:   PetscFunctionReturn(PETSC_SUCCESS);
     97: }
     98: /*@
     99:   PetscDASetOptionsPrefix - Sets the prefix used for searching for all
    100:   PetscDA options in the database.

    102:   Logically Collective

    104:   Input Parameters:
    105: + das - the `PetscDA` context
    106: - p   - the prefix string to prepend to all PetscDA option requests

    108:   Level: advanced

    110: .seealso: `PetscDA`, `PetscDASetFromOptions()`, `PetscDAAppendOptionsPrefix()`, `PetscDAGetOptionsPrefix()`
    111: @*/
    112: PetscErrorCode PetscDASetOptionsPrefix(PetscDA das, const char p[])
    113: {
    114:   PetscFunctionBegin;
    116:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)das, p));
    117:   PetscFunctionReturn(PETSC_SUCCESS);
    118: }

    120: /*@
    121:   PetscDAAppendOptionsPrefix - Appends to the prefix used for searching for all PetscDA options in the database.

    123:   Logically Collective

    125:   Input Parameters:
    126: + das - the `PetscDA` context
    127: - p   - the prefix string to prepend to all `PetscDA` option requests

    129:   Level: advanced

    131: .seealso: `PetscDA`, `PetscDASetFromOptions()`, `PetscDASetOptionsPrefix()`, `PetscDAGetOptionsPrefix()`
    132: @*/
    133: PetscErrorCode PetscDAAppendOptionsPrefix(PetscDA das, const char p[])
    134: {
    135:   PetscFunctionBegin;
    137:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)das, p));
    138:   PetscFunctionReturn(PETSC_SUCCESS);
    139: }

    141: /*@
    142:   PetscDAGetOptionsPrefix - Gets the prefix used for searching for all
    143:   PetscDA options in the database

    145:   Not Collective

    147:   Input Parameter:
    148: . das - the `PetscDA` context

    150:   Output Parameter:
    151: . p - pointer to the prefix string used

    153:   Level: advanced

    155: .seealso: `PetscDA`, `PetscDASetFromOptions()`, `PetscDASetOptionsPrefix()`, `PetscDAAppendOptionsPrefix()`
    156: @*/
    157: PetscErrorCode PetscDAGetOptionsPrefix(PetscDA das, const char *p[])
    158: {
    159:   PetscFunctionBegin;
    161:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)das, p));
    162:   PetscFunctionReturn(PETSC_SUCCESS);
    163: }

    165: /*@
    166:   PetscDACreate - Creates a new `PetscDA` object for data assimilation.

    168:   Collective

    170:   Input Parameter:
    171: . comm - MPI communicator used to create the object

    173:   Output Parameter:
    174: . da_out - newly created `PetscDA` object

    176:   Level: beginner

    178: .seealso: [](ch_da), `PetscDADestroy()`, `PetscDASetType()`, `PetscDASetUp()`
    179: @*/
    180: PetscErrorCode PetscDACreate(MPI_Comm comm, PetscDA *da_out)
    181: {
    182:   PetscDA da;

    184:   PetscFunctionBegin;
    185:   PetscAssertPointer(da_out, 2);
    186:   PetscCall(PetscDAInitializePackage());
    187:   PetscCall(PetscHeaderCreate(da, PETSCDA_CLASSID, "PetscDA", "Data Assimilation", "DA", comm, PetscDADestroy, PetscDAView));
    188:   PetscCall(PetscMemzero(da->ops, sizeof(*da->ops)));
    189:   da->state_size       = 0;
    190:   da->local_state_size = PETSC_DECIDE;
    191:   da->obs_size         = 0;
    192:   da->local_obs_size   = PETSC_DECIDE;
    193:   da->ndof             = 1;
    194:   da->obs_error_var    = NULL;
    195:   da->R                = NULL;
    196:   da->data             = NULL;
    197:   *da_out              = da;
    198:   PetscCall(PetscDASetType(da, PETSCDAETKF));
    199:   PetscFunctionReturn(PETSC_SUCCESS);
    200: }

    202: /*@
    203:   PetscDADestroy - Destroys a `PetscDA` object and releases its resources.

    205:   Collective

    207:   Input Parameter:
    208: . da - pointer to the `PetscDA` object to destroy

    210:   Level: beginner

    212: .seealso: [](ch_da), `PetscDACreate()`
    213: @*/
    214: PetscErrorCode PetscDADestroy(PetscDA *da)
    215: {
    216:   PetscFunctionBegin;
    217:   if (!da || !*da) PetscFunctionReturn(PETSC_SUCCESS);
    219:   if (--((PetscObject)*da)->refct > 0) {
    220:     *da = NULL;
    221:     PetscFunctionReturn(PETSC_SUCCESS);
    222:   }
    223:   PetscTryTypeMethod(*da, destroy);
    224:   PetscCall(PetscHeaderDestroy(da));
    225:   PetscFunctionReturn(PETSC_SUCCESS);
    226: }

    228: /*@
    229:   PetscDASetType - Sets the data assimilation implementation used by a `PetscDA` object.

    231:   Collective

    233:   Input Parameters:
    234: + da   - the `PetscDA` context
    235: - type - name of the implementation (for example `PETSCDAETKF`)

    237:   Level: intermediate

    239: .seealso: [](ch_da), `PetscDAGetType()`, `PetscDARegister()`
    240: @*/
    241: PetscErrorCode PetscDASetType(PetscDA da, PetscDAType type)
    242: {
    243:   PetscErrorCode (*r)(PetscDA);
    244:   PetscBool match;

    246:   PetscFunctionBegin;
    248:   PetscAssertPointer(type, 2);

    250:   PetscCall(PetscObjectTypeCompare((PetscObject)da, type, &match));
    251:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

    253:   PetscCall(PetscDARegisterAll());
    254:   PetscCall(PetscFunctionListFind(PetscDAList, type, &r));
    255:   PetscCheck(r, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDA type: %s", type);

    257:   PetscTryTypeMethod(da, destroy);
    258:   da->ops->destroy = NULL;
    259:   da->data         = NULL;

    261:   PetscCall((*r)(da));
    262:   PetscCall(PetscObjectChangeTypeName((PetscObject)da, type));
    263:   PetscFunctionReturn(PETSC_SUCCESS);
    264: }

    266: /*@
    267:   PetscDAGetType - Gets the name of the implementation currently associated with a `PetscDA`.

    269:   Not Collective

    271:   Input Parameter:
    272: . da - the `PetscDA` context

    274:   Output Parameter:
    275: . type - pointer that will receive the type name (may be `NULL`)

    277:   Level: intermediate

    279: .seealso: [](ch_da), `PetscDASetType()`
    280: @*/
    281: PetscErrorCode PetscDAGetType(PetscDA da, PetscDAType *type)
    282: {
    283:   PetscFunctionBegin;
    285:   if (type) {
    286:     PetscAssertPointer(type, 2);
    287:     *type = ((PetscObject)da)->type_name;
    288:   }
    289:   PetscFunctionReturn(PETSC_SUCCESS);
    290: }

    292: /*@
    293:   PetscDASetFromOptions - Configures a `PetscDA` object from the options database.

    295:   Collective

    297:   Input Parameter:
    298: . da - the `PetscDA` context to set up

    300:   Level: intermediate

    302: .seealso: [](ch_da), `PetscDASetType()`, `PetscObjectOptionsBegin()`
    303: @*/
    304: PetscErrorCode PetscDASetFromOptions(PetscDA da)
    305: {
    306:   char      type_name[256];
    307:   PetscBool type_set;

    309:   PetscFunctionBegin;

    312:   PetscObjectOptionsBegin((PetscObject)da);
    313:   PetscCall(PetscOptionsFList("-petscda_type", "Data assimilation method", "PetscDASetType", PetscDAList, ((PetscObject)da)->type_name, type_name, sizeof(type_name), &type_set));
    314:   if (type_set) PetscCall(PetscDASetType(da, type_name));
    315:   if (da->ops->setfromoptions) PetscCall((*da->ops->setfromoptions)(da, &PetscOptionsObject));
    316:   PetscOptionsEnd();
    317:   PetscFunctionReturn(PETSC_SUCCESS);
    318: }

    320: /*@
    321:   PetscDASetNDOF - Set the number of degrees of freedom per grid point

    323:   Logically Collective

    325:   Input Parameters:
    326: + da   - the `PetscDA` context
    327: - ndof - number of degrees of freedom per grid point (e.g., 2 for shallow water with h and hu)

    329:   Level: intermediate

    331:   Note:
    332:   This must be called before `PetscDASetUp()`. The default is 1 (scalar field).

    334:   Developer Note:
    335:   It is a limitation that each grid point needs the same number of degrees of freedom.

    337: .seealso: `PetscDA`, `PetscDAGetNDOF()`, `PetscDASetUp()`, `PetscDASetSizes()`
    338: @*/
    339: PetscErrorCode PetscDASetNDOF(PetscDA da, PetscInt ndof)
    340: {
    341:   PetscFunctionBegin;
    344:   PetscCheck(ndof > 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "ndof must be positive, got %" PetscInt_FMT, ndof);
    345:   da->ndof = ndof;
    346:   PetscFunctionReturn(PETSC_SUCCESS);
    347: }

    349: /*@
    350:   PetscDAGetNDOF - Get the number of degrees of freedom per grid point

    352:   Not Collective

    354:   Input Parameter:
    355: . da - the `PetscDA` context

    357:   Output Parameter:
    358: . ndof - number of degrees of freedom per grid point

    360:   Level: intermediate

    362: .seealso: `PetscDA`, `PetscDASetNDOF()`
    363: @*/
    364: PetscErrorCode PetscDAGetNDOF(PetscDA da, PetscInt *ndof)
    365: {
    366:   PetscFunctionBegin;
    368:   PetscAssertPointer(ndof, 2);
    369:   *ndof = da->ndof;
    370:   PetscFunctionReturn(PETSC_SUCCESS);
    371: }

    373: /*@
    374:   PetscDASetUp - Allocates internal data structures for a `PetscDA` based on the previously provided sizes.

    376:   Collective

    378:   Input Parameter:
    379: . da - the `PetscDA` context to assemble

    381:   Level: beginner

    383: .seealso: [](ch_da), `PetscDASetSizes()`, `PetscDASetType()`
    384: @*/
    385: PetscErrorCode PetscDASetUp(PetscDA da)
    386: {
    387:   PetscFunctionBegin;
    389:   PetscTryTypeMethod(da, setup);
    390:   PetscFunctionReturn(PETSC_SUCCESS);
    391: }

    393: /*@
    394:   PetscDAView - Views a `PetscDA` and its implementation-specific data structure.

    396:   Collective

    398:   Input Parameters:
    399: + da     - the `PetscDA` context
    400: - viewer - the `PetscViewer` to use (or `NULL` for standard output)

    402:   Level: beginner

    404: .seealso: [](ch_da), `PetscDAViewFromOptions()`
    405: @*/
    406: PetscErrorCode PetscDAView(PetscDA da, PetscViewer viewer)
    407: {
    408:   PetscBool   iascii;
    409:   PetscMPIInt size;

    411:   PetscFunctionBegin;
    413:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)da), &viewer));
    415:   PetscCheckSameComm(da, 1, viewer, 2);

    417:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    418:   if (iascii) {
    419:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)da), &size));
    420:     PetscCall(PetscViewerASCIIPrintf(viewer, "PetscDA Object: %" PetscInt_FMT " MPI process%s\n", (PetscInt)size, size > 1 ? "es" : ""));
    421:     PetscCall(PetscViewerASCIIPrintf(viewer, "  type: %s\n", ((PetscObject)da)->type_name ? ((PetscObject)da)->type_name : "not set"));
    422:     PetscCall(PetscViewerASCIIPrintf(viewer, "  State size: %" PetscInt_FMT "\n", da->state_size));
    423:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Observation size: %" PetscInt_FMT "\n", da->obs_size));
    424:   }

    426:   if (da->ops->view) PetscCall((*da->ops->view)(da, viewer));
    427:   PetscFunctionReturn(PETSC_SUCCESS);
    428: }

    430: /*@
    431:   PetscDAViewFromOptions - Processes command-line options to determine if a `PetscDA` should be viewed.

    433:   Collective

    435:   Input Parameters:
    436: + da   - the `PetscDA` context
    437: . obj  - optional object that provides the prefix for options
    438: - name - option name to check

    440:   Options Database Key:
    441: . -name [viewertype][:...] - option name and values. See `PetscObjectViewFromOptions()` for the possible arguments

    443:   Level: beginner

    445: .seealso: [](ch_da), `PetscDAView()`, `PetscObjectViewFromOptions()`
    446: @*/
    447: PetscErrorCode PetscDAViewFromOptions(PetscDA da, PetscObject obj, const char name[])
    448: {
    449:   PetscFunctionBegin;
    451:   PetscCall(PetscObjectViewFromOptions((PetscObject)da, obj, name));
    452:   PetscFunctionReturn(PETSC_SUCCESS);
    453: }

    455: /*@
    456:   PetscDASetObsErrorVariance - Sets the observation-error variances associated with a `PetscDA`.

    458:   Collective

    460:   Input Parameters:
    461: + da            - the `PetscDA` context
    462: - obs_error_var - vector containing observation error variances (assumes R is a diagonal matrix)

    464:   Notes:
    465:   This function creates or updates both the observation error variance vector and the
    466:   observation error covariance matrix `R`. The matrix `R` is constructed as a diagonal matrix
    467:   with the variances on the diagonal.

    469:   Level: beginner

    471: .seealso: [](ch_da), `PetscDAGetObsErrorVariance()`
    472: @*/
    473: PetscErrorCode PetscDASetObsErrorVariance(PetscDA da, Vec obs_error_var)
    474: {
    475:   PetscFunctionBegin;
    478:   PetscCheck(da->obs_size > 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ORDER, "Call PetscDASetSizes() before PetscDASetObsErrorVariance()");

    480:   /* Create or update observation error variance vector */
    481:   if (!da->obs_error_var) PetscCall(VecDuplicate(obs_error_var, &da->obs_error_var));
    482:   PetscCall(VecCopy(obs_error_var, da->obs_error_var));

    484:   /* Create or update observation error covariance matrix R (p x p) as AIJ matrix
    485:      This is currently initialized as a diagonal matrix, but can be used
    486:      for non-diagonal covariance in the future */
    487:   if (!da->R) {
    488:     PetscCall(MatCreate(PetscObjectComm((PetscObject)da), &da->R));
    489:     PetscCall(MatSetSizes(da->R, da->local_obs_size, da->local_obs_size, da->obs_size, da->obs_size));
    490:     PetscCall(MatSetType(da->R, MATAIJ));
    491:     PetscCall(MatSetFromOptions(da->R));
    492:     PetscCall(MatSetUp(da->R));
    493:   }

    495:   /* Set R as diagonal matrix with variances on diagonal */
    496:   PetscCall(MatZeroEntries(da->R));
    497:   PetscCall(MatDiagonalSet(da->R, da->obs_error_var, INSERT_VALUES));
    498:   PetscCall(MatAssemblyBegin(da->R, MAT_FINAL_ASSEMBLY));
    499:   PetscCall(MatAssemblyEnd(da->R, MAT_FINAL_ASSEMBLY));
    500:   PetscFunctionReturn(PETSC_SUCCESS);
    501: }

    503: /*@
    504:   PetscDAGetObsErrorVariance - Returns a borrowed reference to the observation-error variance vector.

    506:   Not Collective

    508:   Input Parameter:
    509: . da - the `PetscDA` context

    511:   Output Parameter:
    512: . obs_error_var - pointer to the variance vector managed by the `PetscDA`

    514:   Level: beginner

    516: .seealso: [](ch_da), `PetscDASetObsErrorVariance()`
    517: @*/
    518: PetscErrorCode PetscDAGetObsErrorVariance(PetscDA da, Vec *obs_error_var)
    519: {
    520:   PetscFunctionBegin;
    522:   PetscAssertPointer(obs_error_var, 2);
    523:   *obs_error_var = da->obs_error_var;
    524:   PetscFunctionReturn(PETSC_SUCCESS);
    525: }

    527: /*@
    528:   PetscDASetSizes - Sets the state and observation sizes for a `PetscDA`

    530:   Collective

    532:   Input Parameters:
    533: + da         - the `PetscDA` context
    534: . state_size - number of state components
    535: - obs_size   - number of observation components

    537:   Level: beginner

    539:   Developer Note:
    540:   It is not clear this is a good API, shouldn't one provide template vectors for these?

    542: .seealso: [](ch_da), `PetscDAGetSizes()`, `PetscDASetUp()`, `PetscDAEnsembleSetSize()`
    543: @*/
    544: PetscErrorCode PetscDASetSizes(PetscDA da, PetscInt state_size, PetscInt obs_size)
    545: {
    546:   PetscFunctionBegin;

    551:   da->state_size = state_size;
    552:   da->obs_size   = obs_size;
    553:   PetscFunctionReturn(PETSC_SUCCESS);
    554: }

    556: /*@
    557:   PetscDASetLocalSizes - Sets the local state and observation dimensions used by a `PetscDA`.

    559:   Collective

    561:   Input Parameters:
    562: + da               - the `PetscDA` context
    563: . local_state_size - number of local state components (or `PETSC_DECIDE`)
    564: - local_obs_size   - number of local observation components (or `PETSC_DECIDE`)

    566:   Level: beginner

    568: .seealso: [](ch_da), `PetscDASetSizes()`, `PetscDASetUp()`
    569: @*/
    570: PetscErrorCode PetscDASetLocalSizes(PetscDA da, PetscInt local_state_size, PetscInt local_obs_size)
    571: {
    572:   PetscFunctionBegin;
    574:   da->local_state_size = local_state_size;
    575:   da->local_obs_size   = local_obs_size;
    576:   PetscFunctionReturn(PETSC_SUCCESS);
    577: }

    579: /*@
    580:   PetscDAGetSizes - Retrieves the state size and observation size from a `PetscDA`.

    582:   Not Collective

    584:   Input Parameter:
    585: . da - the `PetscDA` context

    587:   Output Parameters:
    588: + state_size - number of state components (may be `NULL`)
    589: - obs_size   - number of observation components (may be `NULL`)

    591:   Level: beginner

    593: .seealso: [](ch_da), `PetscDASetSizes()`
    594: @*/
    595: PetscErrorCode PetscDAGetSizes(PetscDA da, PetscInt *state_size, PetscInt *obs_size)
    596: {
    597:   PetscFunctionBegin;
    599:   if (state_size) PetscAssertPointer(state_size, 2);
    600:   if (obs_size) PetscAssertPointer(obs_size, 3);
    601:   if (state_size) *state_size = da->state_size;
    602:   if (obs_size) *obs_size = da->obs_size;
    603:   PetscFunctionReturn(PETSC_SUCCESS);
    604: }