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_LETKF(PetscDA);
79: /*@
80: PetscDARegisterAll - Registers all data assimilation backends that were compiled in.
82: Not Collective
84: Level: developer
86: .seealso: [](ch_da), `PetscDARegister()`
87: @*/
88: PetscErrorCode PetscDARegisterAll(void)
89: {
90: PetscFunctionBegin;
91: if (PetscDARegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
92: PetscDARegisterAllCalled = PETSC_TRUE;
93: PetscCall(PetscDARegister(PETSCDALETKF, PetscDACreate_LETKF));
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
96: /*@
97: PetscDASetOptionsPrefix - Sets the prefix used for searching for all
98: PetscDA options in the database.
100: Logically Collective
102: Input Parameters:
103: + das - the `PetscDA` context
104: - p - the prefix string to prepend to all PetscDA option requests
106: Level: advanced
108: .seealso: `PetscDA`, `PetscDASetFromOptions()`, `PetscDAAppendOptionsPrefix()`, `PetscDAGetOptionsPrefix()`
109: @*/
110: PetscErrorCode PetscDASetOptionsPrefix(PetscDA das, const char p[])
111: {
112: PetscFunctionBegin;
114: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)das, p));
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: /*@
119: PetscDAAppendOptionsPrefix - Appends to the prefix used for searching for all PetscDA options in the database.
121: Logically Collective
123: Input Parameters:
124: + das - the `PetscDA` context
125: - p - the prefix string to prepend to all `PetscDA` option requests
127: Level: advanced
129: .seealso: `PetscDA`, `PetscDASetFromOptions()`, `PetscDASetOptionsPrefix()`, `PetscDAGetOptionsPrefix()`
130: @*/
131: PetscErrorCode PetscDAAppendOptionsPrefix(PetscDA das, const char p[])
132: {
133: PetscFunctionBegin;
135: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)das, p));
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: /*@
140: PetscDAGetOptionsPrefix - Gets the prefix used for searching for all
141: PetscDA options in the database
143: Not Collective
145: Input Parameter:
146: . das - the `PetscDA` context
148: Output Parameter:
149: . p - pointer to the prefix string used
151: Level: advanced
153: .seealso: `PetscDA`, `PetscDASetFromOptions()`, `PetscDASetOptionsPrefix()`, `PetscDAAppendOptionsPrefix()`
154: @*/
155: PetscErrorCode PetscDAGetOptionsPrefix(PetscDA das, const char *p[])
156: {
157: PetscFunctionBegin;
159: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)das, p));
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: /*@
164: PetscDACreate - Creates a new `PetscDA` object for data assimilation.
166: Collective
168: Input Parameter:
169: . comm - MPI communicator used to create the object
171: Output Parameter:
172: . da_out - newly created `PetscDA` object
174: Level: beginner
176: .seealso: [](ch_da), `PetscDADestroy()`, `PetscDASetType()`, `PetscDASetUp()`
177: @*/
178: PetscErrorCode PetscDACreate(MPI_Comm comm, PetscDA *da_out)
179: {
180: PetscDA da;
182: PetscFunctionBegin;
183: PetscAssertPointer(da_out, 2);
184: PetscCall(PetscDAInitializePackage());
185: PetscCall(PetscHeaderCreate(da, PETSCDA_CLASSID, "PetscDA", "Data Assimilation", "DA", comm, PetscDADestroy, PetscDAView));
186: PetscCall(PetscMemzero(da->ops, sizeof(*da->ops)));
187: da->state_size = 0;
188: da->local_state_size = PETSC_DECIDE;
189: da->obs_size = 0;
190: da->local_obs_size = PETSC_DECIDE;
191: da->ndof = 1;
192: da->obs_error_var = NULL;
193: da->R = NULL;
194: da->data = NULL;
195: *da_out = da;
196: PetscCall(PetscDASetType(da, PETSCDALETKF));
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
200: /*@
201: PetscDADestroy - Destroys a `PetscDA` object and releases its resources.
203: Collective
205: Input Parameter:
206: . da - pointer to the `PetscDA` object to destroy
208: Level: beginner
210: .seealso: [](ch_da), `PetscDACreate()`
211: @*/
212: PetscErrorCode PetscDADestroy(PetscDA *da)
213: {
214: PetscFunctionBegin;
215: if (!da || !*da) PetscFunctionReturn(PETSC_SUCCESS);
217: if (--((PetscObject)*da)->refct > 0) {
218: *da = NULL;
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
221: PetscTryTypeMethod(*da, destroy);
222: PetscCall(PetscHeaderDestroy(da));
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: /*@
227: PetscDASetType - Sets the data assimilation implementation used by a `PetscDA` object.
229: Collective
231: Input Parameters:
232: + da - the `PetscDA` context
233: - type - name of the implementation (for example `PETSCDALETKF`)
235: Level: intermediate
237: .seealso: [](ch_da), `PetscDAGetType()`, `PetscDARegister()`
238: @*/
239: PetscErrorCode PetscDASetType(PetscDA da, PetscDAType type)
240: {
241: PetscErrorCode (*r)(PetscDA);
242: PetscBool match;
244: PetscFunctionBegin;
246: PetscAssertPointer(type, 2);
248: PetscCall(PetscObjectTypeCompare((PetscObject)da, type, &match));
249: if (match) PetscFunctionReturn(PETSC_SUCCESS);
251: PetscCall(PetscDARegisterAll());
252: PetscCall(PetscFunctionListFind(PetscDAList, type, &r));
253: PetscCheck(r, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDA type: %s", type);
255: PetscTryTypeMethod(da, destroy);
256: da->ops->destroy = NULL;
257: da->data = NULL;
259: PetscCall((*r)(da));
260: PetscCall(PetscObjectChangeTypeName((PetscObject)da, type));
261: PetscFunctionReturn(PETSC_SUCCESS);
262: }
264: /*@
265: PetscDAGetType - Gets the name of the implementation currently associated with a `PetscDA`.
267: Not Collective
269: Input Parameter:
270: . da - the `PetscDA` context
272: Output Parameter:
273: . type - pointer that will receive the type name (may be `NULL`)
275: Level: intermediate
277: .seealso: [](ch_da), `PetscDASetType()`
278: @*/
279: PetscErrorCode PetscDAGetType(PetscDA da, PetscDAType *type)
280: {
281: PetscFunctionBegin;
283: if (type) {
284: PetscAssertPointer(type, 2);
285: *type = ((PetscObject)da)->type_name;
286: }
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: /*@
291: PetscDASetFromOptions - Configures a `PetscDA` object from the options database.
293: Collective
295: Input Parameter:
296: . da - the `PetscDA` context to set up
298: Level: intermediate
300: .seealso: [](ch_da), `PetscDASetType()`, `PetscObjectOptionsBegin()`
301: @*/
302: PetscErrorCode PetscDASetFromOptions(PetscDA da)
303: {
304: char type_name[256];
305: PetscBool type_set;
307: PetscFunctionBegin;
310: PetscObjectOptionsBegin((PetscObject)da);
311: PetscCall(PetscOptionsFList("-petscda_type", "Data assimilation method", "PetscDASetType", PetscDAList, ((PetscObject)da)->type_name, type_name, sizeof(type_name), &type_set));
312: if (type_set) PetscCall(PetscDASetType(da, type_name));
313: PetscTryTypeMethod(da, setfromoptions, &PetscOptionsObject);
314: PetscOptionsEnd();
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: /*@
319: PetscDASetNDOF - Set the number of degrees of freedom per grid point
321: Logically Collective
323: Input Parameters:
324: + da - the `PetscDA` context
325: - ndof - number of degrees of freedom per grid point (e.g., 2 for shallow water with h and hu)
327: Level: intermediate
329: Note:
330: This must be called before `PetscDASetUp()`. The default is 1 (scalar field).
332: Developer Note:
333: It is a limitation that each grid point needs the same number of degrees of freedom.
335: .seealso: `PetscDA`, `PetscDAGetNDOF()`, `PetscDASetUp()`, `PetscDASetSizes()`
336: @*/
337: PetscErrorCode PetscDASetNDOF(PetscDA da, PetscInt ndof)
338: {
339: PetscFunctionBegin;
342: PetscCheck(ndof > 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "ndof must be positive, got %" PetscInt_FMT, ndof);
343: da->ndof = ndof;
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: /*@
348: PetscDAGetNDOF - Get the number of degrees of freedom per grid point
350: Not Collective
352: Input Parameter:
353: . da - the `PetscDA` context
355: Output Parameter:
356: . ndof - number of degrees of freedom per grid point
358: Level: intermediate
360: .seealso: `PetscDA`, `PetscDASetNDOF()`
361: @*/
362: PetscErrorCode PetscDAGetNDOF(PetscDA da, PetscInt *ndof)
363: {
364: PetscFunctionBegin;
366: PetscAssertPointer(ndof, 2);
367: *ndof = da->ndof;
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: /*@
372: PetscDASetUp - Allocates internal data structures for a `PetscDA` based on the previously provided sizes.
374: Collective
376: Input Parameter:
377: . da - the `PetscDA` context to assemble
379: Level: beginner
381: .seealso: [](ch_da), `PetscDASetSizes()`, `PetscDASetType()`
382: @*/
383: PetscErrorCode PetscDASetUp(PetscDA da)
384: {
385: PetscFunctionBegin;
387: PetscTryTypeMethod(da, setup);
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: /*@
392: PetscDAView - Views a `PetscDA` and its implementation-specific data structure.
394: Collective
396: Input Parameters:
397: + da - the `PetscDA` context
398: - viewer - the `PetscViewer` to use (or `NULL` for standard output)
400: Level: beginner
402: .seealso: [](ch_da), `PetscDAViewFromOptions()`
403: @*/
404: PetscErrorCode PetscDAView(PetscDA da, PetscViewer viewer)
405: {
406: PetscBool iascii;
407: PetscMPIInt size;
409: PetscFunctionBegin;
411: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)da), &viewer));
413: PetscCheckSameComm(da, 1, viewer, 2);
415: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
416: if (iascii) {
417: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)da), &size));
418: PetscCall(PetscViewerASCIIPrintf(viewer, "PetscDA Object: %" PetscInt_FMT " MPI process%s\n", (PetscInt)size, size > 1 ? "es" : ""));
419: PetscCall(PetscViewerASCIIPrintf(viewer, " type: %s\n", ((PetscObject)da)->type_name ? ((PetscObject)da)->type_name : "not set"));
420: PetscCall(PetscViewerASCIIPrintf(viewer, " State size: %" PetscInt_FMT "\n", da->state_size));
421: PetscCall(PetscViewerASCIIPrintf(viewer, " Observation size: %" PetscInt_FMT "\n", da->obs_size));
422: }
424: PetscTryTypeMethod(da, view, viewer);
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: /*@
429: PetscDAViewFromOptions - Processes command-line options to determine if a `PetscDA` should be viewed.
431: Collective
433: Input Parameters:
434: + da - the `PetscDA` context
435: . obj - optional object that provides the prefix for options
436: - name - option name to check
438: Options Database Key:
439: . -name [viewertype][:...] - option name and values. See `PetscObjectViewFromOptions()` for the possible arguments
441: Level: beginner
443: .seealso: [](ch_da), `PetscDAView()`, `PetscObjectViewFromOptions()`
444: @*/
445: PetscErrorCode PetscDAViewFromOptions(PetscDA da, PetscObject obj, const char name[])
446: {
447: PetscFunctionBegin;
449: PetscCall(PetscObjectViewFromOptions((PetscObject)da, obj, name));
450: PetscFunctionReturn(PETSC_SUCCESS);
451: }
453: /*@
454: PetscDASetObsErrorVariance - Sets the observation-error variances associated with a `PetscDA`.
456: Collective
458: Input Parameters:
459: + da - the `PetscDA` context
460: - obs_error_var - vector containing observation error variances (assumes R is a diagonal matrix)
462: Notes:
463: This function creates or updates both the observation error variance vector and the
464: observation error covariance matrix `R`. The matrix `R` is constructed as a diagonal matrix
465: with the variances on the diagonal.
467: Level: beginner
469: .seealso: [](ch_da), `PetscDAGetObsErrorVariance()`
470: @*/
471: PetscErrorCode PetscDASetObsErrorVariance(PetscDA da, Vec obs_error_var)
472: {
473: PetscFunctionBegin;
476: PetscCheck(da->obs_size > 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ORDER, "Call PetscDASetSizes() before PetscDASetObsErrorVariance()");
478: /* Create or update observation error variance vector */
479: if (!da->obs_error_var) PetscCall(VecDuplicate(obs_error_var, &da->obs_error_var));
480: PetscCall(VecCopy(obs_error_var, da->obs_error_var));
482: /* Create or update observation error covariance matrix R (p x p) as AIJ matrix
483: This is currently initialized as a diagonal matrix, but can be used
484: for non-diagonal covariance in the future */
485: if (!da->R) {
486: PetscCall(MatCreate(PetscObjectComm((PetscObject)da), &da->R));
487: PetscCall(MatSetSizes(da->R, da->local_obs_size, da->local_obs_size, da->obs_size, da->obs_size));
488: PetscCall(MatSetType(da->R, MATAIJ));
489: PetscCall(MatSetFromOptions(da->R));
490: PetscCall(MatSetUp(da->R));
491: }
493: /* Set R as diagonal matrix with variances on diagonal */
494: PetscCall(MatZeroEntries(da->R));
495: PetscCall(MatDiagonalSet(da->R, da->obs_error_var, INSERT_VALUES));
496: PetscCall(MatAssemblyBegin(da->R, MAT_FINAL_ASSEMBLY));
497: PetscCall(MatAssemblyEnd(da->R, MAT_FINAL_ASSEMBLY));
498: PetscFunctionReturn(PETSC_SUCCESS);
499: }
501: /*@
502: PetscDAGetObsErrorVariance - Returns a borrowed reference to the observation-error variance vector.
504: Not Collective
506: Input Parameter:
507: . da - the `PetscDA` context
509: Output Parameter:
510: . obs_error_var - pointer to the variance vector managed by the `PetscDA`
512: Level: beginner
514: .seealso: [](ch_da), `PetscDASetObsErrorVariance()`
515: @*/
516: PetscErrorCode PetscDAGetObsErrorVariance(PetscDA da, Vec *obs_error_var)
517: {
518: PetscFunctionBegin;
520: PetscAssertPointer(obs_error_var, 2);
521: *obs_error_var = da->obs_error_var;
522: PetscFunctionReturn(PETSC_SUCCESS);
523: }
525: /*@
526: PetscDASetSizes - Sets the state and observation sizes for a `PetscDA`
528: Collective
530: Input Parameters:
531: + da - the `PetscDA` context
532: . state_size - number of state components
533: - obs_size - number of observation components
535: Level: beginner
537: Developer Note:
538: It is not clear this is a good API, shouldn't one provide template vectors for these?
540: .seealso: [](ch_da), `PetscDAGetSizes()`, `PetscDASetUp()`, `PetscDAEnsembleSetSize()`
541: @*/
542: PetscErrorCode PetscDASetSizes(PetscDA da, PetscInt state_size, PetscInt obs_size)
543: {
544: PetscFunctionBegin;
549: da->state_size = state_size;
550: da->obs_size = obs_size;
551: PetscFunctionReturn(PETSC_SUCCESS);
552: }
554: /*@
555: PetscDASetLocalSizes - Sets the local state and observation dimensions used by a `PetscDA`.
557: Collective
559: Input Parameters:
560: + da - the `PetscDA` context
561: . local_state_size - number of local state components (or `PETSC_DECIDE`)
562: - local_obs_size - number of local observation components (or `PETSC_DECIDE`)
564: Level: beginner
566: .seealso: [](ch_da), `PetscDASetSizes()`, `PetscDASetUp()`
567: @*/
568: PetscErrorCode PetscDASetLocalSizes(PetscDA da, PetscInt local_state_size, PetscInt local_obs_size)
569: {
570: PetscFunctionBegin;
572: da->local_state_size = local_state_size;
573: da->local_obs_size = local_obs_size;
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: /*@
578: PetscDAGetSizes - Retrieves the state size and observation size from a `PetscDA`.
580: Not Collective
582: Input Parameter:
583: . da - the `PetscDA` context
585: Output Parameters:
586: + state_size - number of state components (may be `NULL`)
587: - obs_size - number of observation components (may be `NULL`)
589: Level: beginner
591: .seealso: [](ch_da), `PetscDASetSizes()`
592: @*/
593: PetscErrorCode PetscDAGetSizes(PetscDA da, PetscInt *state_size, PetscInt *obs_size)
594: {
595: PetscFunctionBegin;
597: if (state_size) PetscAssertPointer(state_size, 2);
598: if (obs_size) PetscAssertPointer(obs_size, 3);
599: if (state_size) *state_size = da->state_size;
600: if (obs_size) *obs_size = da->obs_size;
601: PetscFunctionReturn(PETSC_SUCCESS);
602: }