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: - option - option name to check (may be `NULL`)
440: Level: beginner
442: .seealso: [](ch_da), `PetscDAView()`, `PetscObjectViewFromOptions()`
443: @*/
444: PetscErrorCode PetscDAViewFromOptions(PetscDA da, PetscObject obj, const char option[])
445: {
446: PetscFunctionBegin;
448: PetscCall(PetscObjectViewFromOptions((PetscObject)da, obj, option));
449: PetscFunctionReturn(PETSC_SUCCESS);
450: }
452: /*@
453: PetscDASetObsErrorVariance - Sets the observation-error variances associated with a `PetscDA`.
455: Collective
457: Input Parameters:
458: + da - the `PetscDA` context
459: - obs_error_var - vector containing observation error variances (assumes R is a diagonal matrix)
461: Notes:
462: This function creates or updates both the observation error variance vector and the
463: observation error covariance matrix `R`. The matrix `R` is constructed as a diagonal matrix
464: with the variances on the diagonal.
466: Level: beginner
468: .seealso: [](ch_da), `PetscDAGetObsErrorVariance()`
469: @*/
470: PetscErrorCode PetscDASetObsErrorVariance(PetscDA da, Vec obs_error_var)
471: {
472: PetscFunctionBegin;
475: PetscCheck(da->obs_size > 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ORDER, "Call PetscDASetSizes() before PetscDASetObsErrorVariance()");
477: /* Create or update observation error variance vector */
478: if (!da->obs_error_var) PetscCall(VecDuplicate(obs_error_var, &da->obs_error_var));
479: PetscCall(VecCopy(obs_error_var, da->obs_error_var));
481: /* Create or update observation error covariance matrix R (p x p) as AIJ matrix
482: This is currently initialized as a diagonal matrix, but can be used
483: for non-diagonal covariance in the future */
484: if (!da->R) {
485: PetscCall(MatCreate(PetscObjectComm((PetscObject)da), &da->R));
486: PetscCall(MatSetSizes(da->R, da->local_obs_size, da->local_obs_size, da->obs_size, da->obs_size));
487: PetscCall(MatSetType(da->R, MATAIJ));
488: PetscCall(MatSetFromOptions(da->R));
489: PetscCall(MatSetUp(da->R));
490: }
492: /* Set R as diagonal matrix with variances on diagonal */
493: PetscCall(MatZeroEntries(da->R));
494: PetscCall(MatDiagonalSet(da->R, da->obs_error_var, INSERT_VALUES));
495: PetscCall(MatAssemblyBegin(da->R, MAT_FINAL_ASSEMBLY));
496: PetscCall(MatAssemblyEnd(da->R, MAT_FINAL_ASSEMBLY));
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: /*@
501: PetscDAGetObsErrorVariance - Returns a borrowed reference to the observation-error variance vector.
503: Not Collective
505: Input Parameter:
506: . da - the `PetscDA` context
508: Output Parameter:
509: . obs_error_var - pointer to the variance vector managed by the `PetscDA`
511: Level: beginner
513: .seealso: [](ch_da), `PetscDASetObsErrorVariance()`
514: @*/
515: PetscErrorCode PetscDAGetObsErrorVariance(PetscDA da, Vec *obs_error_var)
516: {
517: PetscFunctionBegin;
519: PetscAssertPointer(obs_error_var, 2);
520: *obs_error_var = da->obs_error_var;
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: /*@
525: PetscDASetSizes - Sets the state and observation sizes for a `PetscDA`
527: Collective
529: Input Parameters:
530: + da - the `PetscDA` context
531: . state_size - number of state components
532: - obs_size - number of observation components
534: Level: beginner
536: Developer Note:
537: It is not clear this is a good API, shouldn't one provide template vectors for these?
539: .seealso: [](ch_da), `PetscDAGetSizes()`, `PetscDASetUp()`, `PetscDAEnsembleSetSize()`
540: @*/
541: PetscErrorCode PetscDASetSizes(PetscDA da, PetscInt state_size, PetscInt obs_size)
542: {
543: PetscFunctionBegin;
548: da->state_size = state_size;
549: da->obs_size = obs_size;
550: PetscFunctionReturn(PETSC_SUCCESS);
551: }
553: /*@
554: PetscDASetLocalSizes - Sets the local state and observation dimensions used by a `PetscDA`.
556: Collective
558: Input Parameters:
559: + da - the `PetscDA` context
560: . local_state_size - number of local state components (or `PETSC_DECIDE`)
561: - local_obs_size - number of local observation components (or `PETSC_DECIDE`)
563: Level: beginner
565: .seealso: [](ch_da), `PetscDASetSizes()`, `PetscDASetUp()`
566: @*/
567: PetscErrorCode PetscDASetLocalSizes(PetscDA da, PetscInt local_state_size, PetscInt local_obs_size)
568: {
569: PetscFunctionBegin;
571: da->local_state_size = local_state_size;
572: da->local_obs_size = local_obs_size;
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: /*@
577: PetscDAGetSizes - Retrieves the state size and observation size from a `PetscDA`.
579: Not Collective
581: Input Parameter:
582: . da - the `PetscDA` context
584: Output Parameters:
585: + state_size - number of state components (may be `NULL`)
586: - obs_size - number of observation components (may be `NULL`)
588: Level: beginner
590: .seealso: [](ch_da), `PetscDASetSizes()`
591: @*/
592: PetscErrorCode PetscDAGetSizes(PetscDA da, PetscInt *state_size, PetscInt *obs_size)
593: {
594: PetscFunctionBegin;
596: if (state_size) PetscAssertPointer(state_size, 2);
597: if (obs_size) PetscAssertPointer(obs_size, 3);
598: if (state_size) *state_size = da->state_size;
599: if (obs_size) *obs_size = da->obs_size;
600: PetscFunctionReturn(PETSC_SUCCESS);
601: }