1: static char help[] = "Tests basic creation and destruction of PetscDA objects, and a simple ETKF analysis step.\n\n";
  2: #include <petscda.h>

  4: int main(int argc, char **argv)
  5: {
  6:   PetscDA     da;
  7:   Mat         H;
  8:   Vec         x_true, y_obs, obs_error_var;
  9:   Vec         x_mean_forecast, x_mean_analysis;
 10:   PetscInt    state_size = 10, obs_size = 10, ensemble_size = 20;
 11:   PetscRandom rng;
 12:   PetscInt    i;
 13:   PetscReal   norm;

 15:   PetscFunctionBeginUser;
 16:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

 18:   /* Create the DA object */
 19:   PetscCall(PetscDACreate(PETSC_COMM_WORLD, &da));
 20:   PetscCall(PetscDASetType(da, PETSCDAETKF));
 21:   PetscCall(PetscDASetSizes(da, state_size, obs_size));
 22:   PetscCall(PetscDAEnsembleSetSize(da, ensemble_size));
 23:   PetscCall(PetscDASetFromOptions(da));
 24:   PetscCall(PetscDASetUp(da));

 26:   /* Initialize random number generator */
 27:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rng));
 28:   PetscCall(PetscRandomSetFromOptions(rng));

 30:   /* Create identity observation matrix H (obs_size x state_size) */
 31:   PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, obs_size, state_size, 1, NULL, 0, NULL, &H));
 32:   PetscCall(MatSetFromOptions(H));
 33:   PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
 34:   PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
 35:   PetscCall(MatShift(H, 1.0));

 37:   /* Create vectors using MatCreateVecs from H */
 38:   PetscCall(MatCreateVecs(H, &x_true, &y_obs));
 39:   PetscCall(VecSet(x_true, 1.0)); /* True state is all 1s */

 41:   PetscCall(VecDuplicate(y_obs, &obs_error_var));
 42:   PetscCall(VecSet(obs_error_var, 0.1)); /* Observation error variance */
 43:   PetscCall(PetscDASetObsErrorVariance(da, obs_error_var));

 45:   /* Create synthetic observation: y = x_true (identity observation, no noise for this simple test) */
 46:   PetscCall(VecCopy(x_true, y_obs));

 48:   /* Initialize ensemble with some spread around 0 (far from truth 1.0) */
 49:   for (i = 0; i < ensemble_size; i++) {
 50:     Vec member;
 51:     PetscCall(VecDuplicate(x_true, &member));
 52:     PetscCall(VecSetRandom(member, rng)); /* Uniform random [0, 1] */
 53:     PetscCall(PetscDAEnsembleSetMember(da, i, member));
 54:     PetscCall(VecDestroy(&member));
 55:   }

 57:   /* Compute forecast mean before analysis */
 58:   PetscCall(VecDuplicate(x_true, &x_mean_forecast));
 59:   PetscCall(PetscDAEnsembleComputeMean(da, x_mean_forecast));

 61:   /* Check forecast error */
 62:   PetscCall(VecAXPY(x_mean_forecast, -1.0, x_true));
 63:   PetscCall(VecNorm(x_mean_forecast, NORM_2, &norm));
 64:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Forecast error norm: %g\n", (double)norm));

 66:   /* Perform Analysis Step */
 67:   PetscCall(PetscDAEnsembleAnalysis(da, y_obs, H));

 69:   /* Compute analysis mean */
 70:   PetscCall(VecDuplicate(x_true, &x_mean_analysis));
 71:   PetscCall(PetscDAEnsembleComputeMean(da, x_mean_analysis));

 73:   /* Check analysis error */
 74:   PetscCall(VecAXPY(x_mean_analysis, -1.0, x_true));
 75:   PetscCall(VecNorm(x_mean_analysis, NORM_2, &norm));
 76:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Analysis error norm: %g\n", (double)norm));

 78:   /* The analysis should move the ensemble closer to the observation (truth) */
 79:   /* Since observation error is small (0.1) and prior spread is ~0.08, it should pull towards observation */

 81:   PetscCall(PetscDAViewFromOptions(da, NULL, "-petscda_view"));

 83:   /* Cleanup */
 84:   PetscCall(MatDestroy(&H));
 85:   PetscCall(VecDestroy(&x_true));
 86:   PetscCall(VecDestroy(&y_obs));
 87:   PetscCall(VecDestroy(&obs_error_var));
 88:   PetscCall(VecDestroy(&x_mean_forecast));
 89:   PetscCall(VecDestroy(&x_mean_analysis));
 90:   PetscCall(PetscRandomDestroy(&rng));
 91:   PetscCall(PetscDADestroy(&da));

 93:   PetscCall(PetscFinalize());
 94:   return 0;
 95: }

 97: /*TEST

 99:   test:
100:     suffix: 1
101:     requires: !complex
102:     args: -petscda_view

104:   test:
105:     suffix: chol
106:     requires: !complex
107:     args: -petscda_view -petscda_ensemble_sqrt_type cholesky

109: TEST*/