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: PetscTryTypeMethod(da, setfromoptions, &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: PetscTryTypeMethod(da, view, 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: }