PetscDA: Data Assimilation

PETSc’s PetscDA object coordinates data assimilation (DA) workflows.

This is new code, please independently verify all results you obtain using it.

Some planned work for PetscDA is available as GitLab Issue #1882

Ensemble-based Data Assimilation

Currently PetscDA supports one ensemble-based assimilator: PETSCDALETKF. It implements the Local Ensemble Transform Kalman Filter, with a none-localization mode that reduces to the classic global ETKF. The framework is extensible to other assimilation techniques.

• LETKF

The PetscDA object owns the ensemble Mat, the observation operator, and the registered forecast/analysis callbacks, so assimilation algorithms work against a single API regardless of how the ensemble is partitioned across MPI ranks or which Mat/Vec backend stores it.

Lifecycle overview

A typical assimilation cycle alternates between forecast propagation and statistical analysis:

1. Initialize a PetscDA context and configure ensemble sizes and data structures.

2. Populate the ensemble state vectors and optional observation-error descriptions.

3. Advance each ensemble member with a model operator supplied by the application.

4. Combine forecasts with observations through PetscDAEnsembleAnalysis() to produce the posterior ensemble.

5. Repeat until the desired simulation horizon is complete, optionally extracting diagnostics after each phase.

Throughout this loop the PetscDA object abstracts the global vectors, scatters, and reductions needed to compute ensemble means, anomalies, and square-root transforms.

Here we introduce a simple example to demonstrate PetscDA usage with the PETSCDALETKF on the Lorenz-96 model. Please read Implementations for more in-depth discussion of the available implementations.

Listing: src/ml/da/tutorials/ex1.c

int main(int argc, char **argv)
{
  Lorenz96Ctx *l95_ctx = NULL, *truth_ctx = NULL;
  DM           da_state;
  PetscDA      da;
  Vec          x0, x_mean, x_forecast;
  Vec          truth_state, rmse_work;
  Vec          observation, obs_noise, obs_error_var;
  Mat          H = NULL; /* Observation operator matrix */
  PetscRandom  rng;
  PetscInt     n = DEFAULT_N, steps = DEFAULT_STEPS, burn = DEFAULT_BURN, obs_freq = DEFAULT_OBS_FREQ;
  PetscInt     random_seed = DEFAULT_RANDOM_SEED, ensemble_size = DEFAULT_ENSEMBLE_SIZE;
  PetscInt     n_stat_steps = 0, obs_count = 0, step, progress_interval;
  PetscReal    F = DEFAULT_F, dt = DEFAULT_DT, obs_error_std = DEFAULT_OBS_ERROR_STD;
  PetscReal    ensemble_init_std = 1; /* Initial ensemble spread */
  PetscReal    rmse_forecast = 0.0, rmse_analysis = 0.0;
  PetscReal    sum_rmse_forecast = 0.0, sum_rmse_analysis = 0.0;

  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &argv, NULL, help));

  /* Parse command-line options */
  PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Lorenz-96 LETKF Quick Example", NULL);
  PetscCall(PetscOptionsInt("-n", "State dimension", "", n, &n, NULL));
  PetscCall(PetscOptionsInt("-steps", "Number of time steps", "", steps, &steps, NULL));
  PetscCall(PetscOptionsInt("-burn", "Burn-in steps excluded from statistics", "", burn, &burn, NULL));
  PetscCall(PetscOptionsInt("-obs_freq", "Observation frequency", "", obs_freq, &obs_freq, NULL));
  PetscCall(PetscOptionsReal("-F", "Forcing parameter", "", F, &F, NULL));
  PetscCall(PetscOptionsReal("-dt", "Time step size", "", dt, &dt, NULL));
  PetscCall(PetscOptionsReal("-obs_error", "Observation error standard deviation", "", obs_error_std, &obs_error_std, NULL));
  PetscCall(PetscOptionsReal("-ensemble_init_std", "Initial ensemble spread standard deviation", "", ensemble_init_std, &ensemble_init_std, NULL));
  PetscCall(PetscOptionsInt("-random_seed", "Random seed for ensemble perturbations", "", random_seed, &random_seed, NULL));
  PetscOptionsEnd();

  /* Validate and constrain parameters */
  PetscCall(ValidateParameters(&n, &steps, &burn, &obs_freq, &ensemble_size, &dt, &F, &obs_error_std));

  /* Calculate progress reporting interval (avoid division by zero) */
  progress_interval = (steps >= PROGRESS_INTERVALS) ? (steps / PROGRESS_INTERVALS) : 1;

  /* Create 1D periodic DM for state space */
  PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, n, 1, 2, NULL, &da_state));
  PetscCall(DMSetFromOptions(da_state));
  PetscCall(DMSetUp(da_state));
  PetscCall(DMDASetUniformCoordinates(da_state, 0.0, (PetscReal)n, 0.0, 0.0, 0.0, 0.0));

  /* Create Lorenz96 context with reusable TS object */
  PetscCall(Lorenz96ContextCreate(da_state, n, F, dt, &l95_ctx));
  PetscCall(Lorenz96ContextCreate(da_state, n, F, dt, &truth_ctx));

  /* Initialize random number generator */
  PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rng));
  PetscCall(PetscRandomSetSeed(rng, (unsigned long)random_seed));
  PetscCall(PetscRandomSetFromOptions(rng));
  PetscCall(PetscRandomSeed(rng));

  /* Initialize state vectors */
  PetscCall(DMCreateGlobalVector(da_state, &x0));
  PetscCall(PetscRandomSetInterval(rng, -.1 * F, .1 * F));
  PetscCall(VecSetRandom(x0, rng));
  PetscCall(PetscRandomSetInterval(rng, 0, 1));

  /* Initialize truth trajectory */
  PetscCall(VecDuplicate(x0, &truth_state));
  PetscCall(VecCopy(x0, truth_state));
  PetscCall(VecDuplicate(x0, &rmse_work));

  /* Spin up truth to get onto attractor */
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Spinning up truth for %" PetscInt_FMT " steps...\n", (PetscInt)SPINUP_STEPS));
  for (PetscInt k = 0; k < SPINUP_STEPS; k++) PetscCall(Lorenz96StepVec(truth_ctx, truth_state));

  /* Initialize observation vectors */
  PetscCall(VecDuplicate(x0, &observation));
  PetscCall(VecDuplicate(x0, &obs_noise));
  PetscCall(VecDuplicate(x0, &obs_error_var));
  PetscCall(VecSet(obs_error_var, obs_error_std * obs_error_std));

  /* Initialize ensemble statistics vectors */
  PetscCall(VecDuplicate(x0, &x_mean));
  PetscCall(VecDuplicate(x0, &x_forecast));

  /* Create identity observation matrix H */
  PetscCall(CreateIdentityObservationMatrix(n, &H));

  /* Create and configure PetscDA for ensemble data assimilation */
  PetscCall(PetscDACreate(PETSC_COMM_WORLD, &da));
  PetscCall(PetscDALETKFSetLocalizationType(da, PETSCDA_LETKF_LOC_NONE));
  PetscCall(PetscDASetSizes(da, n, n));
  PetscCall(PetscDAEnsembleSetSize(da, ensemble_size));
  PetscCall(PetscDASetFromOptions(da));
  PetscCall(PetscDAEnsembleGetSize(da, &ensemble_size));
  PetscCall(PetscDASetUp(da));
  PetscCall(PetscDASetObsErrorVariance(da, obs_error_var));

  /* Initialize ensemble members from spun-up truth state with appropriate spread */
  PetscCall(PetscDAEnsembleInitialize(da, truth_state, ensemble_init_std, rng));

  /* Print configuration summary */
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Lorenz-96 LETKF Example\n"));
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "======================\n"));
  PetscCall(PetscPrintf(PETSC_COMM_WORLD,
                        "  State dimension       : %" PetscInt_FMT "\n"
                        "  Ensemble size         : %" PetscInt_FMT "\n"
                        "  Forcing parameter (F) : %.4f\n"
                        "  Time step (dt)        : %.4f\n"
                        "  Total steps           : %" PetscInt_FMT "\n"
                        "  Burn-in steps         : %" PetscInt_FMT "\n"
                        "  Observation frequency : %" PetscInt_FMT "\n"
                        "  Observation noise std : %.3f\n"
                        "  Ensemble init std     : %.3f\n"
                        "  Random seed           : %" PetscInt_FMT "\n\n",
                        n, ensemble_size, (double)F, (double)dt, steps, burn, obs_freq, (double)obs_error_std, (double)ensemble_init_std, random_seed));

  /* Main assimilation cycle: forecast and analysis steps */
  for (step = 0; step <= steps; step++) {
    PetscReal time = step * dt;

    /* Forecast step: compute ensemble mean and forecast RMSE */
    PetscCall(PetscDAEnsembleComputeMean(da, x_mean));
    PetscCall(VecCopy(x_mean, x_forecast));
    PetscCall(ComputeRMSE(x_forecast, truth_state, rmse_work, n, &rmse_forecast));
    rmse_analysis = rmse_forecast; /* Default to forecast RMSE if no analysis */

    /* Analysis step: assimilate observations when available */
    if (step % obs_freq == 0 && step > 0) {
      /* Generate synthetic noisy observations from truth */
      PetscCall(VecSetRandomGaussian(obs_noise, rng, 0.0, obs_error_std));
      PetscCall(VecWAXPY(observation, 1.0, obs_noise, truth_state));

      /* Perform ETKF analysis with observation matrix H */
      PetscCall(PetscDAEnsembleAnalysis(da, observation, H));

      /* Compute analysis RMSE */
      PetscCall(PetscDAEnsembleComputeMean(da, x_mean));
      PetscCall(ComputeRMSE(x_mean, truth_state, rmse_work, n, &rmse_analysis));
      obs_count++;
    }

    /* Accumulate statistics after burn-in period */
    if (step >= burn) {
      sum_rmse_forecast += rmse_forecast;
      sum_rmse_analysis += rmse_analysis;
      n_stat_steps++;
    }

    /* Progress reporting */
    if ((step % progress_interval == 0) || (step == steps) || (step == 0)) {
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Step %4" PetscInt_FMT ", time %6.3f  RMSE_forecast %.5f  RMSE_analysis %.5f%s\n", step, (double)time, (double)rmse_forecast, (double)rmse_analysis, (step < burn) ? " [burn-in]" : ""));
    }

    /* Propagate ensemble and truth trajectory */
    if (step < steps) {
      PetscCall(PetscDAEnsembleForecast(da, Lorenz96Step, l95_ctx));
      PetscCall(Lorenz96StepVec(truth_ctx, truth_state));
    }
  }

  /* Report final statistics */
  if (n_stat_steps > 0) {
    PetscReal avg_rmse_forecast = sum_rmse_forecast / n_stat_steps;
    PetscReal avg_rmse_analysis = sum_rmse_analysis / n_stat_steps;
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nStatistics (%" PetscInt_FMT " post-burn-in steps):\n", n_stat_steps));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "==================================================\n"));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Mean RMSE (forecast) : %.5f\n", (double)avg_rmse_forecast));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Mean RMSE (analysis) : %.5f\n", (double)avg_rmse_analysis));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Observations used    : %" PetscInt_FMT "\n\n", obs_count));
  } else {
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nWarning: No post-burn-in statistics collected (burn >= steps)\n\n"));
  }

  /* Test VecSetRandomGaussian to verify Gaussian distribution */
  {
    Vec          test_vec;
    PetscInt     test_size = 10000; /* Large sample for statistical testing */
    PetscScalar *array;
    PetscReal    mean_target = 2.0, std_target = 1.5;
    PetscReal    sample_mean = 0.0, sample_variance = 0.0, sample_std;
    PetscReal    skewness = 0.0, kurtosis = 0.0;
    PetscBool    test_gaussian = PETSC_FALSE;

    PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_gaussian", &test_gaussian, NULL));

    if (test_gaussian) {
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n==============================================\n"));
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Testing VecSetRandomGaussian\n"));
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "==============================================\n"));

      /* Create test vector */
      PetscCall(VecCreate(PETSC_COMM_WORLD, &test_vec));
      PetscCall(VecSetSizes(test_vec, PETSC_DECIDE, test_size));
      PetscCall(VecSetFromOptions(test_vec));

      /* Generate Gaussian random numbers */
      PetscCall(VecSetRandomGaussian(test_vec, rng, mean_target, std_target));

      /* Get array for statistical analysis */
      PetscCall(VecGetArray(test_vec, &array));

      /* Compute sample mean */
      for (PetscInt i = 0; i < test_size; i++) sample_mean += PetscRealPart(array[i]);
      sample_mean /= test_size;

      /* Compute sample variance and higher moments */
      for (PetscInt i = 0; i < test_size; i++) {
        PetscReal diff  = PetscRealPart(array[i]) - sample_mean;
        PetscReal diff2 = diff * diff;
        sample_variance += diff2;
        skewness += diff * diff2;
        kurtosis += diff2 * diff2;
      }
      sample_variance /= (test_size - 1);
      sample_std = PetscSqrtReal(sample_variance);

      /* Normalize skewness and kurtosis */
      skewness = (skewness / test_size) / PetscPowReal(sample_std, 3.0);
      kurtosis = (kurtosis / test_size) / PetscPowReal(sample_std, 4.0) - 3.0; /* Excess kurtosis */

      PetscCall(VecRestoreArray(test_vec, &array));

      /* Report results */
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nTarget parameters:\n"));
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Mean      : %.6f\n", (double)mean_target));
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Std Dev   : %.6f\n", (double)std_target));

      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nSample statistics (n=%" PetscInt_FMT "):\n", (PetscInt)test_size));
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Mean      : %.6f (error: %.6f)\n", (double)sample_mean, (double)PetscAbsReal(sample_mean - mean_target)));
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Std Dev   : %.6f (error: %.6f)\n", (double)sample_std, (double)PetscAbsReal(sample_std - std_target)));
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Skewness  : %.6f (expected ~0 for Gaussian)\n", (double)skewness));
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Kurtosis  : %.6f (expected ~0 for Gaussian)\n", (double)kurtosis));

      /* Statistical tests with reasonable tolerances for finite samples */
      PetscReal mean_error     = PetscAbsReal(sample_mean - mean_target);
      PetscReal std_error      = PetscAbsReal(sample_std - std_target);
      PetscReal mean_tolerance = 3.0 * std_target / PetscSqrtReal((PetscReal)test_size); /* 3-sigma rule */
      PetscReal std_tolerance  = 0.1 * std_target;                                       /* 10% tolerance for std dev */
      PetscReal skew_tolerance = 0.1;                                                    /* Skewness should be near 0 */
      PetscReal kurt_tolerance = 0.5;                                                    /* Excess kurtosis should be near 0 */

      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nStatistical tests:\n"));
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Mean test     : %s (error %.6f < tolerance %.6f)\n", mean_error < mean_tolerance ? "PASS" : "FAIL", (double)mean_error, (double)mean_tolerance));
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Std dev test  : %s (error %.6f < tolerance %.6f)\n", std_error < std_tolerance ? "PASS" : "FAIL", (double)std_error, (double)std_tolerance));
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Skewness test : %s (|skewness| %.6f < tolerance %.6f)\n", PetscAbsReal(skewness) < skew_tolerance ? "PASS" : "FAIL", (double)PetscAbsReal(skewness), (double)skew_tolerance));
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Kurtosis test : %s (|kurtosis| %.6f < tolerance %.6f)\n", PetscAbsReal(kurtosis) < kurt_tolerance ? "PASS" : "FAIL", (double)PetscAbsReal(kurtosis), (double)kurt_tolerance));

      /* Overall test result */
      PetscBool all_pass = (mean_error < mean_tolerance) && (std_error < std_tolerance) && (PetscAbsReal(skewness) < skew_tolerance) && (PetscAbsReal(kurtosis) < kurt_tolerance);
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nOverall result: %s\n", all_pass ? "PASS - Distribution is Gaussian" : "FAIL - Distribution may not be Gaussian"));
      PetscCall(PetscPrintf(PETSC_COMM_WORLD, "==============================================\n\n"));

      PetscCall(VecDestroy(&test_vec));
    }
  }

  PetscCall(PetscDAView(da, PETSC_VIEWER_STDOUT_WORLD));

  /* Cleanup */
  PetscCall(MatDestroy(&H));
  PetscCall(VecDestroy(&x_forecast));
  PetscCall(VecDestroy(&x_mean));
  PetscCall(VecDestroy(&obs_error_var));
  PetscCall(VecDestroy(&obs_noise));
  PetscCall(VecDestroy(&observation));
  PetscCall(VecDestroy(&rmse_work));
  PetscCall(VecDestroy(&truth_state));
  PetscCall(VecDestroy(&x0));
  PetscCall(PetscDADestroy(&da));
  PetscCall(DMDestroy(&da_state));
  PetscCall(Lorenz96ContextDestroy(&l95_ctx));
  PetscCall(Lorenz96ContextDestroy(&truth_ctx));
  PetscCall(PetscRandomDestroy(&rng));

  PetscCall(PetscFinalize());
return 0;}

Creating a PetscDA context

Create, configure, and destroy a PetscDA object with the standard PETSc object lifecycle:

PetscDA da;
PetscCall(PetscDACreate(PETSC_COMM_WORLD, &da));
PetscCall(PetscDASetType(da, PETSCDALETKF));
PetscCall(PetscDASetSizes(da, state_size, obs_size));
PetscCall(PetscDAEnsembleSetSize(da, ensemble_size));
PetscCall(PetscDASetFromOptions(da));
PetscCall(PetscDASetUp(da));

/* ... data assimilation loop ... */

PetscCall(PetscDADestroy(&da));

To create a PetscDA instance, call PetscDACreate():

PetscDACreate(MPI_Comm comm, PetscDA *da);

To choose an implementation type, call

PetscDASetType(PetscDA da, PetscDAType type);

or use the command-line option -petscda_type name; details regarding the available implementations are presented in Implementations.

PetscDASetSizes() records the global state dimension and the number of observations:

PetscDASetSizes(PetscDA da, PetscInt state_size, PetscInt obs_size)

For MPI-parallel runs where the state or observation vectors are distributed, the local partition sizes can be set explicitly with

PetscDASetLocalSizes(PetscDA da, PetscInt local_state_size, PetscInt local_obs_size);

Pass PETSC_DECIDE for either argument to let PETSc choose the partition automatically.

PetscDAEnsembleSetSize() records the number of ensemble members requested:

PetscDAEnsembleSetSize(PetscDA da, PetscInt ensemble_size)

After having set these options, call PetscDASetFromOptions() to apply any command-line overrides, then PetscDASetUp() to allocate internal storage:

PetscDASetFromOptions(PetscDA da);
PetscDASetUp(PetscDA da);

Finally, after the user is done using the DA object, destroy it with

PetscDADestroy(PetscDA *da);

Degrees of freedom per grid point

When the state vector is laid out on a structured grid, the number of physical degrees of freedom per grid point can be recorded with

PetscDASetNDOF(PetscDA da, PetscInt ndof);
PetscDAGetNDOF(PetscDA da, PetscInt *ndof);

This metadata is used by some implementations and viewers; the default is 1.

Managing ensembles

PetscDA stores ensemble members as PETSc Vec objects and exposes convenience helpers to access them safely.

Individual ensemble members can be read or overwritten using a get/restore ownership pattern analogous to MatDenseGetColumnVec:

/* Read-only view of member i; must be paired with Restore */
PetscDAEnsembleGetMember(PetscDA da, PetscInt i, Vec *member);
PetscDAEnsembleRestoreMember(PetscDA da, PetscInt i, Vec *member);

/* Inject an externally created vector into slot i */
PetscDAEnsembleSetMember(PetscDA da, PetscInt i, Vec member);

PetscDAEnsembleGetMember() / PetscDAEnsembleRestoreMember() map ensemble indices to Vec handles that participate in PETSc’s reference counting. PetscDAEnsembleSetMember() lets applications inject externally created vectors into specific slots, which is useful when importing state snapshots from disk or another solver component.

Ensemble statistics are computed with:

/* Form the sample mean across all members */
PetscDAEnsembleComputeMean(PetscDA da, Vec mean);

/* Return a tall-and-skinny Mat whose columns are the mean-subtracted,
   normalized anomalies X = (E - mean) / sqrt(m-1) */
PetscDAEnsembleComputeAnomalies(PetscDA da, Vec mean, Mat *anomalies);

PetscDAEnsembleComputeAnomalies() returns a dense Mat whose columns are the mean-subtracted ensemble anomalies. Many square-root filters use this matrix to construct low-rank covariance factorizations. Pass NULL for mean to have the function compute it internally.

Observation error

Observation-error variances (or more general descriptions) are supplied through PetscDASetObsErrorVariance(). The associated vector is assumed to follow the global observation ordering; PetscDAGetObsErrorVariance() returns the stored object for later inspection or reuse.

PetscDASetObsErrorVariance(PetscDA da, Vec obs_error_var);
PetscDAGetObsErrorVariance(PetscDA da, Vec *obs_error_var);

Analysis step

PetscDAEnsembleAnalysis() performs the assimilation step given an observation vector and a linear observation operator matrix:

/* observation - Vec of length P (number of observations) */
/* H           - Mat of size P x N (observation operator mapping state to observation space) */
PetscCall(PetscDAEnsembleAnalysis(da, observation, H));

H is a Mat (typically sparse AIJ) that maps the N-dimensional state vector to the P-dimensional observation space: y ≈ H*x. PetscDAAnalysis() handles all ensemble reductions, gain computations, and posterior updates.

Forecast step

PetscDAEnsembleForecast() wraps the forecast step. The user supplies a function that advances the entire ensemble matrix in one call:

Calling sequence for model:

• ensemble - the ensemble matrix whose columns are the states to be advanced in place

• ctx - the context for the model function

A model that can advance the whole ensemble at once (e.g. a vectorized RHS, a Kokkos-resident propagator, or a single time integrator over a stacked state) writes directly to the columns of ensemble. A model that can only advance one state at a time, such as a TS-driven step, loops over the columns itself. Use MatDenseGetColumnVec() / MatDenseRestoreColumnVec() (the read-write pair) for the common case where the model reads the current column as an initial condition and writes the advanced state back; reach for MatDenseGetColumnVecRead() / MatDenseGetColumnVecWrite() only when the kernel is truly read-only or write-only (e.g. a write-only resampler):

/* Prototype for model forecast M(X) */
PetscErrorCode ModelForecast(Mat ensemble, PetscCtx ctx) {
  PetscInt n, j;

  PetscFunctionBeginUser;
  PetscCall(MatGetSize(ensemble, NULL, &n));
  for (j = 0; j < n; j++) {
    Vec col;

    PetscCall(MatDenseGetColumnVec(ensemble, j, &col));
    /* Advance col by dt (e.g., step a TS object) */
    PetscCall(MatDenseRestoreColumnVec(ensemble, j, &col));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscCall(PetscDAEnsembleForecast(da, ModelForecast, ctx));

No MatAssemblyBegin()/MatAssemblyEnd() call is needed after MatDenseRestoreColumnVec(): the dense Get/Restore pair is itself the assembled write path, and the matrix stays assembled across it.

ModelForecast can call PETSc time integrators (TS: Scalable ODE and DAE Solvers), nonlinear solvers (SNES: Nonlinear Solvers), or bespoke device kernels.

Inflation

Covariance inflation counteracts ensemble collapse by artificially widening the prior spread before the analysis step. The inflation factor \(\rho \ge 1\) scales each anomaly so that the effective prior covariance becomes \(\rho^2 P^f\). Set and retrieve the factor with

PetscDAEnsembleSetInflation(PetscDA da, PetscReal inflation);
PetscDAEnsembleGetInflation(PetscDA da, PetscReal *inflation);

or at runtime with -petscda_ensemble_inflation value (default: 1.0, i.e. no inflation).

Implementations

The only built-in PetscDA implementation is the Local Ensemble Transform Kalman Filter (PETSCDALETKF, selected with -petscda_type letkf). Custom types can be registered with PetscDARegister() and selected at runtime with -petscda_type name.

LETKF

The Local ETKF (PETSCDALETKF, -petscda_type letkf) is the default implementation. It performs the analysis update locally around each grid point, enabling scalable assimilation on large domains by avoiding the global ensemble covariance matrix. With -petscda_letkf_localization_type none it reduces to the classic global ETKF (Algorithm 6.4 in [ABN16]), a deterministic square-root update that avoids stochastic perturbations.

The reduced-space T-matrix is factored via a symmetric eigendecomposition \(T = V D V^T\), so the square root used in the ensemble transform is the symmetric \(T^{-1/2} = V D^{-1/2} V^T\). The symmetric form minimizes the rotation of the prior ensemble (preserving member continuity across analysis cycles) and is the only square root that is consistent across overlapping local domains under localization.

LETKF-specific configuration:

/* Distance-based localization: pick a kernel, set the radius, supply
   per-dimension state coordinates and the observation operator H.
   The localization matrix Q is built lazily on the first analysis. */
PetscDALETKFSetLocalizationType(PetscDA da, PetscDALETKFLocalizationType type);
PetscDALETKFSetLocalizationRadius(PetscDA da, PetscReal radius);
PetscDALETKFSetLocalizationCoordinates(PetscDA da, const Vec xyz[3], const PetscReal bd[3], Mat H);

The built-in -petscda_letkf_localization_type values gaspari_cohn, gaussian, and boxcar are available in every PETSc build; the none type disables localization and is mathematically equivalent to global ETKF. The default kernel is gaspari_cohn. The localization matrix Q is built with the same matrix type as the internal observation-error covariance matrix R (assembled from the variance vector passed to PetscDASetObsErrorVariance()): a Kokkos backend is used when R has type MATAIJKOKKOS, otherwise a CPU analysis path is used. Select the kernel at runtime with -petscda_letkf_localization_type (none|gaspari_cohn|gaussian|boxcar).

Both the CPU and Kokkos analysis paths run multi-rank. The unlocalized fast path (PETSCDA_LETKF_LOC_NONE) replicates the m × m T matrix on PETSC_COMM_SELF on every rank: each rank computes its slice of the gram \(S^T S\) and the projection \(S^T \delta\) with local BLAS, the contributions are combined with MPIU_Allreduce(), and the symmetric eigendecomposition runs identically on every rank. The localized per-vertex path uses a backend-agnostic observation scatter that gathers the unique observations referenced by each rank’s rows of Q into a sequential work vector before the per-vertex updates.

Command-line options

The PetscDA object obeys standard PETSc options parsing. Commonly used switches include:

• -petscda_type name – select a registered PetscDA implementation (letkf).

• -petscda_ensemble_inflation value – set the covariance inflation factor (default: 1.0).

• -petscda_view – inspect ensemble metadata and internal sizes.

Because PetscDA participates in the PETSc object registry, any prefix applied with PetscDASetOptionsPrefix() scopes these options.

Viewing and monitoring

PetscDAView() and PetscDAViewFromOptions() expose ensemble sizing, observation dimensions, and implementation-specific diagnostics . Views can be directed to ASCII, HDF5, or custom PetscViewer targets, enabling lightweight instrumentation of assimilation experiments. For advanced profiling, the PetscDA package registers with PETSc’s logging infrastructure via PetscDAInitializePackage()/PetscDAFinalizePackage(), so standard -log_view outputs will include assimilation breakdown.

Further reading

• TS: Scalable ODE and DAE Solvers discusses PETSc time integrators that can supply the forecast operator passed to PetscDAEnsembleForecast().

• Vectors and Parallel Data documents vector assembly and parallel data management for the state and observation spaces.

• SNES: Nonlinear Solvers outlines nonlinear solvers that often participate in observation or model operators.

• DM Basics provides background on distributed mesh infrastructure that can coexist with PetscDA-managed ensembles.

[ABN16]

M. Asch, M. Bocquet, and M. Nodet. Data Assimilation: Methods, Algorithms, and Applications. SIAM, 2016. doi:10.1137/1.9781611974546.