1: #pragma once

  3: #include <petscmat.h>

  5: /* MANSEC = ML */
  6: /* SUBMANSEC = PetscDA */

  8: /*S
  9:    PetscDA - Abstract PETSc object that manages data assimilation

 11:    Level: beginner

 13:    Notes:
 14:    This is new code, please independently verify all results you obtain using it.
 15:    Planned work for `PetscDA` is tracked on GitLab as Issue #1882.
 16:    The only ensemble-based assimilator currently supplied is `PETSCDALETKF`.

 18: .seealso: [](ch_da), `PetscDAType`, `PETSCDALETKF`, `PetscDACreate()`, `PetscDASetType()`,
 19:           `PetscDASetSizes()`, `PetscDAEnsembleSetSize()`, `PetscDAEnsembleAnalysis()`, `PetscDAEnsembleForecast()`,
 20:           `PetscDADestroy()`, `PetscDAView()`
 21: S*/
 22: typedef struct _p_PetscDA *PetscDA;

 24: /*E
 25:   PetscDALETKFLocalizationType - Type of localization kernel used by `PETSCDALETKF`

 27:   Values:
 28: + `PETSCDA_LETKF_LOC_NONE`         - No localization. Each vertex sees every observation with weight one;
 29:                                      the per-vertex loop reduces to a single global analysis (the classic ETKF).
 30: . `PETSCDA_LETKF_LOC_GASPARI_COHN` - Gaspari-Cohn fifth-order piecewise rational kernel with compact support at twice the radius
 31: . `PETSCDA_LETKF_LOC_GAUSSIAN`     - Gaussian kernel exp(-d^2 / (2 r^2)) truncated at twice the radius
 32: - `PETSCDA_LETKF_LOC_BOXCAR`       - Uniform weight one inside the radius, zero outside

 34:   Options Database Keys:
 35: . -petscda_letkf_localization_type (none|gaspari_cohn|gaussian|boxcar) - select the localization kernel at run time

 37:   Level: intermediate

 39: .seealso: [](ch_da), `PETSCDALETKF`, `PetscDALETKFSetLocalizationType()`, `PetscDALETKFGetLocalizationType()`,
 40:           `PetscDALETKFSetLocalizationRadius()`, `PetscDALETKFSetLocalizationCoordinates()`
 41: E*/
 42: typedef enum {
 43:   PETSCDA_LETKF_LOC_NONE,
 44:   PETSCDA_LETKF_LOC_GASPARI_COHN,
 45:   PETSCDA_LETKF_LOC_GAUSSIAN,
 46:   PETSCDA_LETKF_LOC_BOXCAR,
 47:   PETSCDA_LETKF_LOC_NUM_TYPES
 48: } PetscDALETKFLocalizationType;

 50: /*J
 51:   PetscDAType - String with the name of a PETSc data assimilation method

 53:   Level: beginner

 55: .seealso: [](ch_da), `PetscDA`, `PetscDASetType()`, `PETSCDALETKF`
 56: J*/
 57: typedef const char *PetscDAType;
 58: #define PETSCDALETKF "letkf"

 60: PETSC_EXTERN PetscErrorCode PetscDAInitializePackage(void);
 61: PETSC_EXTERN PetscErrorCode PetscDAFinalizePackage(void);

 63: PETSC_EXTERN PetscErrorCode PetscDARegister(const char[], PetscErrorCode (*)(PetscDA));
 64: PETSC_EXTERN PetscErrorCode PetscDARegisterAll(void);

 66: PETSC_EXTERN PetscErrorCode PetscDACreate(MPI_Comm, PetscDA *);
 67: PETSC_EXTERN PetscErrorCode PetscDADestroy(PetscDA *);
 68: PETSC_EXTERN PetscErrorCode PetscDASetType(PetscDA, PetscDAType);
 69: PETSC_EXTERN PetscErrorCode PetscDAGetType(PetscDA, PetscDAType *);
 70: PETSC_EXTERN PetscErrorCode PetscDAView(PetscDA, PetscViewer);
 71: PETSC_EXTERN PetscErrorCode PetscDAViewFromOptions(PetscDA, PetscObject, const char[]);
 72: PETSC_EXTERN PetscErrorCode PetscDASetFromOptions(PetscDA);

 74: PETSC_EXTERN PetscErrorCode PetscDASetSizes(PetscDA, PetscInt, PetscInt);
 75: PETSC_EXTERN PetscErrorCode PetscDASetLocalSizes(PetscDA, PetscInt, PetscInt);
 76: PETSC_EXTERN PetscErrorCode PetscDAGetSizes(PetscDA, PetscInt *, PetscInt *);
 77: PETSC_EXTERN PetscErrorCode PetscDASetNDOF(PetscDA, PetscInt);
 78: PETSC_EXTERN PetscErrorCode PetscDAGetNDOF(PetscDA, PetscInt *);
 79: PETSC_EXTERN PetscErrorCode PetscDASetUp(PetscDA);

 81: PETSC_EXTERN PetscErrorCode PetscDASetObsErrorVariance(PetscDA, Vec);
 82: PETSC_EXTERN PetscErrorCode PetscDAGetObsErrorVariance(PetscDA, Vec *);

 84: PETSC_EXTERN PetscErrorCode PetscDASetOptionsPrefix(PetscDA, const char[]);
 85: PETSC_EXTERN PetscErrorCode PetscDAAppendOptionsPrefix(PetscDA, const char[]);
 86: PETSC_EXTERN PetscErrorCode PetscDAGetOptionsPrefix(PetscDA, const char *[]);

 88: PETSC_EXTERN PetscErrorCode PetscDAEnsembleSetSize(PetscDA, PetscInt);
 89: PETSC_EXTERN PetscErrorCode PetscDAEnsembleGetSize(PetscDA, PetscInt *);
 90: PETSC_EXTERN PetscErrorCode PetscDAEnsembleSetInflation(PetscDA, PetscReal);
 91: PETSC_EXTERN PetscErrorCode PetscDAEnsembleGetInflation(PetscDA, PetscReal *);

 93: PETSC_EXTERN PetscErrorCode PetscDAEnsembleGetMember(PetscDA, PetscInt, Vec *);
 94: PETSC_EXTERN PetscErrorCode PetscDAEnsembleRestoreMember(PetscDA, PetscInt, Vec *);
 95: PETSC_EXTERN PetscErrorCode PetscDAEnsembleSetMember(PetscDA, PetscInt, Vec);

 97: PETSC_EXTERN PetscErrorCode PetscDAEnsembleComputeMean(PetscDA, Vec);
 98: PETSC_EXTERN PetscErrorCode PetscDAEnsembleComputeAnomalies(PetscDA, Vec, Mat *);
 99: PETSC_EXTERN PetscErrorCode PetscDAEnsembleAnalysis(PetscDA, Vec, Mat);

101: /*S
102:   PetscDAEnsembleForecastFn - A prototype of an ensemble forecast model callback passed to `PetscDAEnsembleForecast()`

104:   Calling Sequence:
105: + ensemble - ensemble matrix advanced in place; on entry its columns are the current members, on
106:              return its columns are the members at the next time
107: - ctx      - [optional] user-defined model context

109:   Level: intermediate

111:   Note:
112:   Per-member integrators that read and overwrite each column (the common case, e.g. a `TS`-driven step)
113:   iterate with `MatDenseGetColumnVec()` / `MatDenseRestoreColumnVec()`. Use the read-only or write-only
114:   variants (`MatDenseGetColumnVecRead()` / `MatDenseGetColumnVecWrite()`) only when the kernel truly
115:   does not need the other half (e.g. a write-only resampler).

117: .seealso: [](ch_da), `PetscDA`, `PETSCDALETKF`, `PetscDAEnsembleForecast()`
118: S*/
119: PETSC_EXTERN_TYPEDEF typedef PetscErrorCode PetscDAEnsembleForecastFn(Mat ensemble, PetscCtx ctx);

121: PETSC_EXTERN PetscErrorCode PetscDAEnsembleForecast(PetscDA, PetscDAEnsembleForecastFn *, PetscCtx);
122: PETSC_EXTERN PetscErrorCode PetscDAEnsembleInitialize(PetscDA, Vec, PetscReal, PetscRandom);

124: PETSC_EXTERN PetscErrorCode PetscDAEnsembleComputeNormalizedInnovationMatrix(Mat, Vec, Vec, PetscInt, PetscScalar, Mat);
125: PETSC_EXTERN PetscErrorCode PetscDAEnsembleTFactor(PetscDA, Mat);
126: PETSC_EXTERN PetscErrorCode PetscDAEnsembleApplyTInverse(PetscDA, Vec, Vec);
127: PETSC_EXTERN PetscErrorCode PetscDAEnsembleApplySqrtTInverse(PetscDA, Mat, Mat);

129: PETSC_EXTERN PetscErrorCode PetscDALETKFSetLocalizationRadius(PetscDA, PetscReal);
130: PETSC_EXTERN PetscErrorCode PetscDALETKFGetLocalizationRadius(PetscDA, PetscReal *);
131: PETSC_EXTERN PetscErrorCode PetscDALETKFSetLocalizationType(PetscDA, PetscDALETKFLocalizationType);
132: PETSC_EXTERN PetscErrorCode PetscDALETKFGetLocalizationType(PetscDA, PetscDALETKFLocalizationType *);
133: PETSC_EXTERN PetscErrorCode PetscDALETKFSetLocalizationCoordinates(PetscDA, const Vec[3], const PetscReal[3], Mat);
134: PETSC_EXTERN PetscErrorCode PetscDALETKFResetLocalization(PetscDA);