1: static char help[] = "Tests basic creation and destruction of PetscDA objects, and a simple LETKF (NONE-localization) 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:   PetscReal   norm;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

 80:   PetscCall(PetscDAView(da, PETSC_VIEWER_STDOUT_WORLD));

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

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

 96: /*TEST

 98:   test:
 99:     suffix: 1
100:     requires: !complex

102: TEST*/