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: /*@
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: /*@
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: /*@
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, pass `NULL` to use the options prefix of `da`
436: - name - option name to check
438: Options Database Key:
439: . -name viewer_specification - See `PetscOptionsCreateViewer()` for the values of `viewer_specification`
441: Level: beginner
443: Note:
444: This checks the options database, creates the viewer on-the-fly, uses it and then destroys it. Hence it should not be called in heavily used routines,
445: rather `PetscOptionsCreateViewer()` should be used to construct the viewer once which can then be utilized in the heavily used routine.
447: .seealso: [](ch_da), `PetscDAView()`, `PetscObjectViewFromOptions()`, `PetscOptionsCreateViewer()`
448: @*/
449: PetscErrorCode PetscDAViewFromOptions(PetscDA da, PetscObject obj, const char name[])
450: {
451: PetscFunctionBegin;
453: PetscCall(PetscObjectViewFromOptions((PetscObject)da, obj, name));
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: /*@
458: PetscDASetObsErrorVariance - Sets the observation-error variances associated with a `PetscDA`.
460: Collective
462: Input Parameters:
463: + da - the `PetscDA` context
464: - obs_error_var - vector containing observation error variances (assumes R is a diagonal matrix)
466: Notes:
467: This function creates or updates both the observation error variance vector and the
468: observation error covariance matrix `R`. The matrix `R` is constructed as a diagonal matrix
469: with the variances on the diagonal.
471: Level: beginner
473: .seealso: [](ch_da), `PetscDAGetObsErrorVariance()`
474: @*/
475: PetscErrorCode PetscDASetObsErrorVariance(PetscDA da, Vec obs_error_var)
476: {
477: PetscFunctionBegin;
480: PetscCheck(da->obs_size > 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ORDER, "Call PetscDASetSizes() before PetscDASetObsErrorVariance()");
482: /* Create or update observation error variance vector */
483: if (!da->obs_error_var) PetscCall(VecDuplicate(obs_error_var, &da->obs_error_var));
484: PetscCall(VecCopy(obs_error_var, da->obs_error_var));
486: /* Create or update observation error covariance matrix R (p x p) as AIJ matrix
487: This is currently initialized as a diagonal matrix, but can be used
488: for non-diagonal covariance in the future */
489: if (!da->R) {
490: PetscCall(MatCreate(PetscObjectComm((PetscObject)da), &da->R));
491: PetscCall(MatSetSizes(da->R, da->local_obs_size, da->local_obs_size, da->obs_size, da->obs_size));
492: PetscCall(MatSetType(da->R, MATAIJ));
493: PetscCall(MatSetFromOptions(da->R));
494: PetscCall(MatSetUp(da->R));
495: }
497: /* Set R as diagonal matrix with variances on diagonal */
498: PetscCall(MatZeroEntries(da->R));
499: PetscCall(MatDiagonalSet(da->R, da->obs_error_var, INSERT_VALUES));
500: PetscCall(MatAssemblyBegin(da->R, MAT_FINAL_ASSEMBLY));
501: PetscCall(MatAssemblyEnd(da->R, MAT_FINAL_ASSEMBLY));
502: PetscFunctionReturn(PETSC_SUCCESS);
503: }
505: /*@
506: PetscDAGetObsErrorVariance - Returns a borrowed reference to the observation-error variance vector.
508: Not Collective
510: Input Parameter:
511: . da - the `PetscDA` context
513: Output Parameter:
514: . obs_error_var - pointer to the variance vector managed by the `PetscDA`
516: Level: beginner
518: .seealso: [](ch_da), `PetscDASetObsErrorVariance()`
519: @*/
520: PetscErrorCode PetscDAGetObsErrorVariance(PetscDA da, Vec *obs_error_var)
521: {
522: PetscFunctionBegin;
524: PetscAssertPointer(obs_error_var, 2);
525: *obs_error_var = da->obs_error_var;
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
529: /*@
530: PetscDASetSizes - Sets the state and observation sizes for a `PetscDA`
532: Collective
534: Input Parameters:
535: + da - the `PetscDA` context
536: . state_size - number of state components
537: - obs_size - number of observation components
539: Level: beginner
541: Developer Note:
542: It is not clear this is a good API, shouldn't one provide template vectors for these?
544: .seealso: [](ch_da), `PetscDAGetSizes()`, `PetscDASetUp()`, `PetscDAEnsembleSetSize()`
545: @*/
546: PetscErrorCode PetscDASetSizes(PetscDA da, PetscInt state_size, PetscInt obs_size)
547: {
548: PetscFunctionBegin;
553: da->state_size = state_size;
554: da->obs_size = obs_size;
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }
558: /*@
559: PetscDASetLocalSizes - Sets the local state and observation dimensions used by a `PetscDA`.
561: Collective
563: Input Parameters:
564: + da - the `PetscDA` context
565: . local_state_size - number of local state components (or `PETSC_DECIDE`)
566: - local_obs_size - number of local observation components (or `PETSC_DECIDE`)
568: Level: beginner
570: .seealso: [](ch_da), `PetscDASetSizes()`, `PetscDASetUp()`
571: @*/
572: PetscErrorCode PetscDASetLocalSizes(PetscDA da, PetscInt local_state_size, PetscInt local_obs_size)
573: {
574: PetscFunctionBegin;
576: da->local_state_size = local_state_size;
577: da->local_obs_size = local_obs_size;
578: PetscFunctionReturn(PETSC_SUCCESS);
579: }
581: /*@
582: PetscDAGetSizes - Retrieves the state size and observation size from a `PetscDA`.
584: Not Collective
586: Input Parameter:
587: . da - the `PetscDA` context
589: Output Parameters:
590: + state_size - number of state components (may be `NULL`)
591: - obs_size - number of observation components (may be `NULL`)
593: Level: beginner
595: .seealso: [](ch_da), `PetscDASetSizes()`
596: @*/
597: PetscErrorCode PetscDAGetSizes(PetscDA da, PetscInt *state_size, PetscInt *obs_size)
598: {
599: PetscFunctionBegin;
601: if (state_size) PetscAssertPointer(state_size, 2);
602: if (obs_size) PetscAssertPointer(obs_size, 3);
603: if (state_size) *state_size = da->state_size;
604: if (obs_size) *obs_size = da->obs_size;
605: PetscFunctionReturn(PETSC_SUCCESS);
606: }