This program was contributed by Jonas Plank <[email protected]>.
It comes without any warranty or support by its authors or the authors of deal.II.

This program is part of the deal.II code gallery and consists of the following files (click to inspect):

Pictures from this code gallery program

Annotated version of README.md

Multilevel Monte Carlo for random Darcy flow

This code is an implementation of the Multilevel Monte Carlo (MLMC) approach described in the amazing paper by Cliffe et al. [3]. While many Darcy flow solvers in porous media research rely on the Finite Volume Method (FVM), this project specifically utilizes the Finite Element Method (FEM), leveraging the deal.II library for spatial discretization.

Motivation

Numerical simulations rely on input parameters such as material data, boundary conditions, and geometry. These parameters are almost always derived from measurements and are therefore subject to uncertainty. Some PDEs, such as the Kardar-Parisi-Zhang equation [1], even incorporate uncertainty directly into the governing equation. This uncertainty propagates through the simulation; consequently, quantifying its impact on the solution is of great interest.

However, this can be a costly undertaking. The standard Monte Carlo estimator, which simply solves the equation \(N\) times with different realizations of the random noise and then averages the \(N\) solutions, converges at a rate of

\[ \mathcal{O}(N^{-\frac{1}{2}}) \]

, potentially requiring tens of thousands of samples for a reasonable estimate of the mean. This is often infeasible for complex problems where computing a single sample on a fine mesh may take hours.

This is where the Multilevel Monte Carlo (MLMC) estimator is utilized. The method exploits a simple identity: to estimate the expectation of a problem on a fine mesh

\[ \mathcal{E}[Q_M]\]

, we can use the linearity of expectation to expand the estimator as a telescoping sum:

\[ \mathcal{E}[Q_M] = \mathcal{E}[Q_0] + \sum_{l=1}^{M} \mathcal{E}[Q_l - Q_{l-1}] \]

.

One run on level

\[ M-1 \]

is significantly cheaper than one run on level

\[ M \]

. Furthermore, the difference between

\[ Q_M \]

and

\[ Q_{M-1} \]

is likely to be small, as the solutions will be similar when using the same realizations of the random field. By running the discretized PDE on several mesh levels, we can allocate more samples to the inexpensive coarse levels and fewer samples to the expensive fine levels. This significantly reduces the total computational effort.

A prominent application for this method is Darcy's Law [2], which describes flow through porous media. In this context, randomness is introduced via the hydraulic conductivity tensor.

Problem statement

Let the infinite probability space be

\[(\Omega, \mathcal{F}, \mathbb{P})\]

and the domain be

\[D \subset \mathbb{R}^d, d=1,2,3 \]

. We model the hydraulic conductivity as a random field

\[ k = k(x, \omega)\]

on

\[D \times \Omega \]

with a defined mean and covariance. The governing equation is:

\[ -\nabla \cdot ( k(x,\omega)\nabla p(x, \omega)) = f(x). \]

We assume the right-hand side

\[f(x) \]

is deterministic. Solving this general form is challenging; therefore, we assume

\[ k(x, \omega) \]

follows a log-normal distribution. This involves replacing the conductivity tensor with a scalar-valued field whose logarithm is Gaussian, which guarantees that

\[ k > 0 \]

almost surely [3]. Additionally, for an exponential kernel, analytical solutions for the integral eigenvalue problem are available.

We approximate this random field using a Karhunen-Loève Expansion and set a pressure difference of 1 between the left and right boundaries.

It is easy to observe that the actual PDE code is very similar to the one demonstrated in step-5. This is on purpose. In recent years, invasive techniques have fallen out of favor since they, as the name suggests, require reprogramming of existing finite element codes.

Results

The following plots illustrate the convergence rate of our MLMC implementation. Convergence Plots

In the left image, the variance significantly decreases with each additional level. This supports the observation that solutions become increasingly similar as the discretization is refined for the same samples. The central plot shows the change in the mean with each level. The fact that the mean continues to shift substantially at levels 2 and 3 suggests that the initial mesh may have been too coarse. This is further supported by the right plot, which shows a sharp drop in the total number of samples required to reach the target tolerance. While the number of samples required for levels 1 and 2 is similar—again suggesting an overly coarse initial discretization—there is a significant reduction in samples for the higher levels. We conclude that the implementation successfully achieves the expected MLMC convergence behavior.

These plots can be generated using the Python script located in utils/.

Finally, a single realization of the solution is shown below. One realization of the solution

References

[1] https://en.wikipedia.org/wiki/Kardar-Parisi-Zhang_equation [2] https://en.wikipedia.org/wiki/Darcy%27s_law [3] Cliffe, K.A., Giles, M.B., Scheichl, R. et al. Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients. Comput. Visual Sci. 14, 3 (2011). https://doi.org/10.1007/s00791-011-0160-x

Annotated version of include/KL_expansion.h

  /* -----------------------------------------------------------------------------
   *
   * SPDX-License-Identifier: LGPL-2.1-or-later
   * Copyright (C) 2026 by Jonas Plank
   *
   * This file is part of the deal.II code gallery.
   *
   * -----------------------------------------------------------------------------
   */
  #pragma once
  #include <deal.II/base/point.h>
  #include <deal.II/base/function.h>
  
  #include <iostream>
  #include <vector>
  
  namespace RandomField
  {
      using namespace dealii;
      /* We assume here that our auto-covariance function is exponential, but not Gaussian.
      This is both modeling choice and mathematically reasonable in this case.
      We also assume a variance of 1 here.  */
      template <int dim>
      class KLExpansion
      {
          public:
          KLExpansion(const unsigned int n_terms, double L, double l, double mu); 
  
          double compute_kl_expansion(const Point<dim>& p, std::vector<double> &samples);
          
          private:
  

computes our coefficients

computes our eigenvalues based on the previously computed frequencies

void compute_lambda_i();

computes norm values such that our eigenfunctions are normed to one.

void compute_alpha_i();

compute the frequencies

void compute_omega_i();

functions we need for the computation of omega_i

          double f_even(double sol);     
          double f_odd(double sol);       
          double grad_f_even(double sol); 
          double grad_f_odd(double sol);  
  

newton since equation is nonlinear

          void newton_even(unsigned int index);  
          void newton_odd(unsigned int index);   
  

number of KL expansion terms

unsigned int n_terms_;

domain length

double domain_length_;

correlation length

double correlation_length_;

important parameters for the construction of the KL-Expansion in x-direction

          std::vector<double> omega_i;
          std::vector<double> lambda_i;
          std::vector<double> alpha_i;
  

constant mean of our field.

          double mu_;
      };
  }

Annotated version of include/mlmc.h

  /* -----------------------------------------------------------------------------
   *
   * SPDX-License-Identifier: LGPL-2.1-or-later
   * Copyright (C) 2026 by Jonas Plank
   *
   * This file is part of the deal.II code gallery.
   *
   * -----------------------------------------------------------------------------
   */
  #pragma once
  #include <vector>
  #include <random>
  
  namespace MultilevelMonteCarlo
  {
      template <int dim>
      class MLMC 
      {
          public:
          MLMC(unsigned int oneD_samples);                                   
  

generates the required number of samples

std::vector<double> generate_samples();

store the samples for the computation of mean and variance

          void add_sample(double rvalue);                                     
  
          double compute_mean();
  
          double compute_variance();

removes all samples after reached convergence on a level

          void clear_samples();
  
          private:
          std::mt19937 rng;                       // Mersenne Twister engine
          std::normal_distribution<double> dist;  // Normal distribution
          unsigned int num_samples;
  

parameters for error computation

          std::vector<double> results;
      };
  }
  

Annotated version of include/random_darcy.h

  /* -----------------------------------------------------------------------------
   *
   * SPDX-License-Identifier: LGPL-2.1-or-later
   * Copyright (C) 2026 by Jonas Plank
   *
   * This file is part of the deal.II code gallery.
   *
   * -----------------------------------------------------------------------------
   */
  #pragma once
  #include <deal.II/base/quadrature_lib.h>
  #include <deal.II/dofs/dof_handler.h>
  #include <deal.II/dofs/dof_tools.h>
  #include <deal.II/fe/fe_values.h>
  #include <deal.II/grid/tria.h>
  #include <deal.II/grid/grid_generator.h>
  #include <deal.II/lac/dynamic_sparsity_pattern.h>
  #include <deal.II/lac/full_matrix.h>
  #include <deal.II/lac/sparse_matrix.h>
  #include <deal.II/lac/vector.h>
  #include <deal.II/numerics/data_out.h>
  #include <deal.II/numerics/vector_tools.h>
  #include <deal.II/fe/fe_q.h>
  #include <deal.II/grid/grid_out.h>
  #include <deal.II/lac/affine_constraints.h>
  #include <deal.II/grid/grid_refinement.h>
  #include <deal.II/numerics/error_estimator.h>
  #include <deal.II/base/function.h>
  #include <fstream>
  
  #include "random_permeability.h"
  
  namespace Discretization
  {
      using namespace dealii;
  
      
      template <int dim>
      class RandomDarcy
      {
      public:
          RandomDarcy();                                                                   
  

Allows switching between coarse and fine triangulations with virtually no code duplication.

          void set_tria(bool fine);                                                       
  
          void generate_mesh(double domain_length);

Implementation follows the logic of deal.II tutorial step-5.

          void setup_system();                                                            
  
          void assemble_system(RandomField::RandomPermeability<dim>& random_constant); 
  
          void solve();

If we are on the first level, we only want to refine the fine triangulation.

void refine_grid(bool firstRun);

On the first run, we label the output as "coarse" since both meshes are identical. After the first level, the meshes diverge, and we output the finer mesh.

void output_results(bool firstRun, unsigned int level);

level

unsigned int level

Definition grid_out.cc:4642

This is our Quantity of Interest (QoI). It is defined as: \(K_{eff} = - \int_{\Gamma_{right}} k \frac{\partial p}{\partial x_1} dx_2\)

          double compute_Keff(RandomField::RandomPermeability<dim>& permeability);
  
          private:

define two triangulations

          Triangulation<dim> coarse_tria;
          Triangulation<dim> fine_tria;
  
          const FE_Q<dim> fe;
          DoFHandler<dim> dof_handler;
  
          AffineConstraints<double> constraints;
  
          SparseMatrix<double> system_matrix;
          SparsityPattern      sparsity_pattern;
          
          Vector<double> solution;
          Vector<double> system_rhs;
      };
  }

Annotated version of include/random_permeability.h

  /* -----------------------------------------------------------------------------
   *
   * SPDX-License-Identifier: LGPL-2.1-or-later
   * Copyright (C) 2026 by Jonas Plank
   *
   * This file is part of the deal.II code gallery.
   *
   * -----------------------------------------------------------------------------
   */
  #pragma once
  #include "KL_expansion.h"
  
  namespace RandomField
  {
      using namespace dealii;
  
      /*
      This class is really our interface between the KL expansion 
      and the FEM code. */
      template <int dim>
      class RandomPermeability : public Function<dim>
      {
          public:
          RandomPermeability(std::vector<double> &first_sample, unsigned int n_terms, double domain_length, double correlation_length, double mu);  
  

overwrite the samples

void overwrite_samples(const std::vector<double> &next_sample);

our evaluator

          double value(const Point<dim>& p);  
  
          private:
          KLExpansion<dim> kl_expansion;
          std::vector<double> samples;
      };
  }
  

Annotated version of main.cc

  #include "include/random_permeability.h"
  #include "include/random_darcy.h"
  #include "include/mlmc.h"
  
  #include <fstream>
  #include <iostream>
  
  int main()
  {
      using namespace dealii;
  
      std::ofstream outfile("mlmc_results.txt");
      if (!outfile.is_open()) {
          std::cerr << "Error: Could not open mlmc_results.txt for writing!" << std::endl;
          return 1;
      }
  

Run parameters

      double domain_length = 10.0;
  
      /* Describes how strongly different points are correlated. 
      A higher correlation length means the random field is more "spread out." 
      A shorter correlation length increases local variance, which requires more KL terms 
      to capture the higher frequencies. This is similar to a Fourier series, where a 
      function with high-frequency oscillations requires more terms for a sufficient approximation. */ 
      double correlation_length = 5.0;
  

The number of terms chosen for 1D. For higher dimensions, we use oneD_samples^dim.

unsigned int oneD_samples = 20;

The number of levels was not chosen arbitrarily. Convergence is usually so fast that using more than 4–5 levels is rarely necessary.

unsigned int levels = 5;

We used a constant mean here. While this was a specific choice for this case, in practice, most mean functions are spatially varying rather than constant.

      double mu = 1.0;
  
      /* This tolerance should not be confused with the tolerances seen in standard numerical 
      convergence studies or linear solvers. Those low tolerances (on the order of 1e-8) 
      would be realistically unachievable here; instead, we use a more moderate tolerance 
      that keeps the error "small enough." A typical value used in research is 1e-2 to 1e-4 
      at most. A larger value was chosen here for demonstration purposes, but the results 
      are still quite acceptable. Note that this value is squared later, so the actual 
      tolerance used to abort the runs is smaller. */
      double tolerance = 1e-1; 
  

In general, these values should be determined by a pilot run. I chose these specific values to keep the implementation simple.

      std::vector<unsigned int> runs_per_level{20000, 10000, 5000, 2500, 1000, 250};
  
      MultilevelMonteCarlo::MLMC<2> mlmc(oneD_samples);
      std::vector<double> first_sample{};
  
      /* The construction of this class also computes the KL expansion. 
      We compute the expansion once and then reuse it by swapping out the samples. */
      RandomField::RandomPermeability<2> permeability(first_sample, oneD_samples, domain_length, correlation_length, mu);
  
      Discretization::RandomDarcy<2> random_darcy;
  
      random_darcy.generate_mesh(domain_length);
  
      double global_mean = 0.0;
  
  
  
      for(unsigned int i = 0; i<=levels; i++)
      {
          std::cout << "Level:" << i << std::endl;
          for(unsigned int j = 0; j<runs_per_level[i]; j++)
          {
              permeability.overwrite_samples(mlmc.generate_samples());
              

coarse run

              random_darcy.set_tria(false);
              random_darcy.setup_system();
              random_darcy.assemble_system(permeability);
              random_darcy.solve();
              double Keff_coarse = random_darcy.compute_Keff(permeability);
              if( i == 0)
              {
                  mlmc.add_sample(Keff_coarse);
              }
              
  

fine run

              if(i!=0)
              {
                  random_darcy.set_tria(true);
                  random_darcy.setup_system();
                  random_darcy.assemble_system(permeability);
                  random_darcy.solve();
                  double Keff_fine = random_darcy.compute_Keff(permeability);
  
                  mlmc.add_sample(Keff_fine-Keff_coarse);
  
              }
  
              if( j >=2)
              {
                  double var = mlmc.compute_variance();
                  
                  if (var / (j+1) < tolerance * tolerance)
                  {
                      global_mean+= mlmc.compute_mean();
                      std::cout << "number of samples for level" << i << ":" << j << std::endl;
                      random_darcy.output_results(i==0, i);
  
                      outfile << "FINAL_STATS Level:" << i << " Mean:" << mlmc.compute_mean() 
                              << " Var:" << mlmc.compute_variance() << " Samples:" << j << std::endl;
                      break;
                  }
              }
              
  
          }
          random_darcy.refine_grid(i==0);
          mlmc.clear_samples();
          std::cout <<"global mean" << global_mean << std::endl;
      }
  
      outfile.close();
  
      
      return 0;
  }
  

Annotated version of source/KL_expansion.cc

  /* -----------------------------------------------------------------------------
   *
   * SPDX-License-Identifier: LGPL-2.1-or-later
   * Copyright (C) 2026 by Jonas Plank
   *
   * This file is part of the deal.II code gallery.
   *
   * -----------------------------------------------------------------------------
   */
  #include "../include/KL_expansion.h"
  #include <cmath>
  #include <math.h>
  
  template <int dim>
  RandomField::KLExpansion<dim>::KLExpansion(const unsigned int n_terms, double domain_length, double correlation_length, double mu)
  : n_terms_(n_terms)
  , domain_length_(domain_length)
  , correlation_length_(correlation_length)
  , mu_(mu)
  {
      compute_omega_i();
      compute_alpha_i();
      compute_lambda_i();
  }
  
  template <int dim>
  double RandomField::KLExpansion<dim>::compute_kl_expansion(const Point<dim>& p, std::vector<double> &samples)
  {
      double val = mu_;
      if constexpr (dim == 1)
      {
          for(unsigned int i = 0; i<n_terms_; i++)
          {

even

              if((i+1)%2 == 0)
              {
                  val += std::sqrt(lambda_i[i])*samples[i]*alpha_i[i]*std::sin(omega_i[i]*(p[0]-domain_length_/2));
              }

odd

              else 
              {
                  val += std::sqrt(lambda_i[i])*samples[i]*alpha_i[i]*std::cos(omega_i[i]*(p[0]-domain_length_/2));
              }
          }
      }
      if constexpr (dim == 2)
      {
          for(unsigned int i = 0; i<n_terms_; i++)
          {
              for(unsigned int j = 0; j<n_terms_; j++)
              {

both even

                  if((i+1)%2==0 && (j+1)%2==0)
                  {
                      val+=std::sqrt(lambda_i[i]*lambda_i[j])*samples[i*n_terms_+j]*alpha_i[i]*alpha_i[j]*std::sin(omega_i[i]*(p[0]-domain_length_/2))*std::sin(omega_i[j]*(p[1]-domain_length_/2));
                  }

both odd

                  else if((i+1)%2==1 && (j+1)%2==1)
                  {
                      val+=std::sqrt(lambda_i[i]*lambda_i[j])*samples[i*n_terms_+j]*alpha_i[i]*alpha_i[j]*std::cos(omega_i[i]*(p[0]-domain_length_/2))*std::cos(omega_i[j]*(p[1]-domain_length_/2));
                  }

i odd

                  else if((i+1)%2==1 && (j+1)%2==0)
                  {
                      val+=std::sqrt(lambda_i[i]*lambda_i[j])*samples[i*n_terms_+j]*alpha_i[i]*alpha_i[j]*std::cos(omega_i[i]*(p[0]-domain_length_/2))*std::sin(omega_i[j]*(p[1]-domain_length_/2));
                  }

j odd

                  else if((i+1)%2==0 && (j+1)%2==1)
                  {
                      val+=std::sqrt(lambda_i[i]*lambda_i[j])*samples[i*n_terms_+j]*alpha_i[i]*alpha_i[j]*std::sin(omega_i[i]*(p[0]-domain_length_/2))*std::cos(omega_i[j]*(p[1]-domain_length_/2));
                  }
  
              }
          }
      }
  
      return val;
  }
  
  template <int dim>
  void RandomField::KLExpansion<dim>::compute_omega_i()
  {
      for(unsigned int i = 1; i<=n_terms_; i++)
      {
          if(i%2 == 0)
          {
              newton_even(i);
          }
          else
          {
              newton_odd(i);
          }
      }
  }
  
  template <int dim>
  void RandomField::KLExpansion<dim>::compute_alpha_i()
  {
      for(unsigned int i = 1; i<=n_terms_; i++)
      {
          if(i%2 == 0)
          {
              double w_i = omega_i[i-1];
              double sqrt_val = domain_length_/2 -std::sin(w_i*domain_length_)/(2*w_i);
              alpha_i.push_back(1/(std::sqrt(sqrt_val)));
          }
          else
          {
              double w_i = omega_i[i-1];
              double sqrt_val = domain_length_/2 +std::sin(w_i*domain_length_)/(2*w_i);
              alpha_i.push_back(1/(std::sqrt(sqrt_val)));
          }
      }
  }
  
  template <int dim>
  void RandomField::KLExpansion<dim>::compute_lambda_i()
  {
      for(unsigned int i = 1; i<=n_terms_; i++)
      {
          double w_i = omega_i[i-1];
          double lambda = 2*correlation_length_/(1+w_i*w_i*correlation_length_*correlation_length_);
          lambda_i.push_back(lambda);
      }
  
  }
  
  template <int dim>
  void RandomField::KLExpansion<dim>::newton_even(unsigned int index)
  {
      double a = M_PI / domain_length_ * (index - 1);
      double b = M_PI / domain_length_ * index;
  
      double x = a + 0.1 * (b - a);  
  
      for (int i = 0; i < 50; ++i)
      {
          double fx  = f_even(x);
          double dfx = grad_f_even(x);
  
          if (std::abs(dfx) < 1e-14) break; 
  
          double step = fx / dfx;
          double x_new = x - step;
  
          int backtrack_count = 0;
          while ((x_new <= a || x_new >= b) && backtrack_count < 10) 
          {
              step *= 0.5;
              x_new = x - step;
              backtrack_count++;
          }
  
          if (x_new <= a || x_new >= b) break;
  
  
          if (!std::isfinite(x_new)) break;
          if (std::abs(x_new - x) < 1e-12) 
          {
              x = x_new; 
              break;
          }
  
          x = x_new;
      }
      omega_i.push_back(x);
  }
  
  template <int dim>
  void RandomField::KLExpansion<dim>::newton_odd(unsigned int index)
  {
      double x = M_PI / domain_length_ * (static_cast<double>(index) - 0.2);
  
      for (int i = 0; i < 50; ++i)
      {
          double fx  = f_odd(x);
          double dfx = grad_f_odd(x);
  
          if (std::abs(dfx) < 1e-12) break;
  
          double x_new = x - fx / dfx;
  
          if (std::abs(x_new - x) < 1e-10) break;
          if (!std::isfinite(x_new)) break;
  
          x = x_new;
      }
      omega_i.push_back(x);
  }
  
  
  template <int dim>
  double RandomField::KLExpansion<dim>::f_odd(double x)
  {
      return 1.0 / correlation_length_ - x * std::tan(x * domain_length_ / 2.0);
  }
  
  template <int dim>
  double RandomField::KLExpansion<dim>::f_even(double x)
  {
      return (1.0 / correlation_length_) * std::tan(x * domain_length_ / 2.0) + x;
  }
  
  template <int dim>
  double RandomField::KLExpansion<dim>::grad_f_odd(double x)
  {
      const double L = domain_length_;
      const double t = std::tan(x * L / 2.0);
      const double sec2 = 1.0 / std::cos(x * L / 2.0);
      return -t - x * (L / 2.0) * sec2 * sec2;
  }
  
  template <int dim>
  double RandomField::KLExpansion<dim>::grad_f_even(double x)
  {
      const double L = domain_length_;
      const double sec2 = 1.0 / std::cos(x * L / 2.0);
      return (1.0 / correlation_length_) * (L / 2.0) * sec2 * sec2 + 1.0;
  }
  
  template class RandomField::KLExpansion<1>;
  template class RandomField::KLExpansion<2>;
  
  

Annotated version of source/mlmc.cc

  /* -----------------------------------------------------------------------------
   *
   * SPDX-License-Identifier: LGPL-2.1-or-later
   * Copyright (C) 2026 by Jonas Plank
   *
   * This file is part of the deal.II code gallery.
   *
   * -----------------------------------------------------------------------------
   */
  
  #include "../include/mlmc.h"
  #include <cmath>
  
  template <int dim>
  MultilevelMonteCarlo::MLMC<dim>::MLMC(unsigned int oneD_samples)
  : rng(std::random_device{}())   
  , dist(0.0, 1.0)                 // mean = 0, stddev = 1
  {
      if constexpr (dim == 1) num_samples = oneD_samples;
      else if constexpr (dim == 2) num_samples = oneD_samples*oneD_samples;
  }
  
  template <int dim>
  std::vector<double> MultilevelMonteCarlo::MLMC<dim>::generate_samples()
  {
      std::vector<double> samples(num_samples);
  
      for (auto &s : samples)
      {
          s = dist(rng);
      }
      return samples;
  }
  
  template<int dim>
  void MultilevelMonteCarlo::MLMC<dim>::add_sample(double rvalue)
  {
      results.push_back(rvalue);
  }
  
  template <int dim>
  double MultilevelMonteCarlo::MLMC<dim>::compute_mean()
  {
      double mean = 0.0;
      for(unsigned int i = 0; i<results.size(); i++)
      {
          mean+=results[i];
      }
      return mean/results.size();
  }
  
  template <int dim>
  double MultilevelMonteCarlo::MLMC<dim>::compute_variance()
  {
      double mean = compute_mean();
      double var = 0.0;
      for(unsigned int i = 0; i<results.size(); i++)
      {
          var += std::pow((results[i]-mean),2);
      }
  
      return var / (results.size() - 1);
  }
  
  template <int dim>
  void MultilevelMonteCarlo::MLMC<dim>::clear_samples()
  {
      results.clear();
  }
  
  template class MultilevelMonteCarlo::MLMC<1>;
  template class MultilevelMonteCarlo::MLMC<2>;
  

Annotated version of source/random_darcy.cc

  /* -----------------------------------------------------------------------------
   *
   * SPDX-License-Identifier: LGPL-2.1-or-later
   * Copyright (C) 2026 by Jonas Plank
   *
   * This file is part of the deal.II code gallery.
   *
   * -----------------------------------------------------------------------------
   */
  
  #include "../include/random_darcy.h"
  #include <deal.II/lac/sparse_direct.h>
  
  template <int dim>
  Discretization::RandomDarcy<dim>::RandomDarcy()
  : fe(1)
  , dof_handler(coarse_tria)
  {}
  
  template <int dim>
  void Discretization::RandomDarcy<dim>::set_tria(bool fine)
  {
      if(fine)
      {
          dof_handler.reinit(fine_tria);
      }
      else 
      {
          dof_handler.reinit(coarse_tria);
      }
  }
  
  template <int dim>
  void Discretization::RandomDarcy<dim>::generate_mesh(double domain_length)
  {
      GridGenerator::hyper_cube(coarse_tria, 0, domain_length, true);
      GridGenerator::hyper_cube(fine_tria, 0, domain_length, true);
  
      coarse_tria.refine_global(3);
      fine_tria.refine_global(3);      
  }
  
  template <int dim>
  void Discretization::RandomDarcy<dim>::setup_system()
  {
      dof_handler.distribute_dofs(fe);
   
      solution.reinit(dof_handler.n_dofs());
      system_rhs.reinit(dof_handler.n_dofs());
      
      constraints.clear();
      DoFTools::make_hanging_node_constraints(dof_handler, constraints);
  

Left boundary (indicator 0): u = 1.0

      VectorTools::interpolate_boundary_values(dof_handler,
                                               0,
                                               Functions::ConstantFunction<dim>(1.0),
                                               constraints);
  

Right boundary (indicator 1): u = 0.0

      VectorTools::interpolate_boundary_values(dof_handler,
                                               1,
                                               Functions::ConstantFunction<dim>(0.0),
                                               constraints);
      
      constraints.close();
      
      DynamicSparsityPattern dsp(dof_handler.n_dofs());
      DoFTools::make_sparsity_pattern(dof_handler,
                                      dsp,
                                      constraints,
                                      /*keep_constrained_dofs = */ false);
      
      sparsity_pattern.copy_from(dsp);
      
      system_matrix.reinit(sparsity_pattern);
  }
  
  template <int dim>
  void Discretization::RandomDarcy<dim>::assemble_system(RandomField::RandomPermeability<dim>& permeability)
  {
      const QGauss<dim> quadrature_formula(fe.degree + 1);
   
    FEValues<dim> fe_values(fe,
                            quadrature_formula,
                            update_values | update_gradients |
                              update_quadrature_points | update_JxW_values);
   
    const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
   
    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
    Vector<double>     cell_rhs(dofs_per_cell);
   
    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
   
    for (const auto &cell : dof_handler.active_cell_iterators())
      {
        fe_values.reinit(cell);
   
        cell_matrix = 0;
        cell_rhs    = 0;
   
        for (const unsigned int q_index : fe_values.quadrature_point_indices())
          {
            const double current_coefficient =
              permeability.value(fe_values.quadrature_point(q_index));
            for (const unsigned int i : fe_values.dof_indices())
              {
                for (const unsigned int j : fe_values.dof_indices())
                  cell_matrix(i, j) +=
                    (current_coefficient *              // a(x_q)
                     fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
                     fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
                     fe_values.JxW(q_index));           // dx
   
                cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
                                1.0 *                               // f(x)
                                fe_values.JxW(q_index));            // dx
              }
          }
   
        cell->get_dof_indices(local_dof_indices);
        constraints.distribute_local_to_global(
          cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
      }
      
  }
  
  template <int dim>
  void Discretization::RandomDarcy<dim>::solve()
  {
      SparseDirectUMFPACK A_direct;
   
      solution = system_rhs;
      A_direct.solve(system_matrix, solution);
  
      constraints.distribute(solution);
  }
  
  template <int dim>
  void Discretization::RandomDarcy<dim>::refine_grid(bool firstRun)
  {
      fine_tria.refine_global(1);
      if(!firstRun)
      {
          coarse_tria.refine_global(1);
      }
  }
  
  template <int dim>
  void Discretization::RandomDarcy<dim>::output_results(bool firstRun, unsigned int level)
  {
      DataOut<dim> data_out;
    
      data_out.attach_dof_handler(dof_handler);
      data_out.add_data_vector(solution, "solution");
    
      data_out.build_patches();
    
      std::ofstream output(firstRun ? "solutionCoarse" + std::to_string(level) + ".vtk" : "solutionFine" + std::to_string(level) + ".vtk");
      data_out.write_vtk(output);
  }
  
  template <int dim>
  double Discretization::RandomDarcy<dim>::compute_Keff(RandomField::RandomPermeability<dim>& permeability)
  {
    const QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
    FEFaceValues<dim> fe_face_values(fe,
                                         face_quadrature_formula,
                                         update_gradients | update_normal_vectors |
                                     update_JxW_values | update_quadrature_points);
  
    std::vector<Tensor<1, dim>> solution_gradients(face_quadrature_formula.size());
    double Keff = 0.0;
  
    for (const auto &cell : dof_handler.active_cell_iterators())
      { for (const auto face_no : cell->face_indices())
        {
          if (cell->face(face_no)->at_boundary() &&
                (cell->face(face_no)->boundary_id() == 1))
              {
                fe_face_values.reinit(cell, face_no);
                fe_face_values.get_function_gradients(solution, solution_gradients);
                for (const unsigned int q_index : fe_face_values.quadrature_point_indices())
                  {
                    const double current_coefficient = permeability.value(fe_face_values.quadrature_point(q_index));
                    Keff -= current_coefficient*solution_gradients[q_index][0]*fe_face_values.JxW(q_index); 
                  }
              }
  
        }
      }
  
      return Keff;
  }
  
  template class Discretization::RandomDarcy<1>;
  template class Discretization::RandomDarcy<2>;
  

Annotated version of source/random_permeability.cc

  /* -----------------------------------------------------------------------------
   *
   * SPDX-License-Identifier: LGPL-2.1-or-later
   * Copyright (C) 2026 by Jonas Plank
   *
   * This file is part of the deal.II code gallery.
   *
   * -----------------------------------------------------------------------------
   */
  #include "../include/random_permeability.h"
  
  template <int dim>
  RandomField::RandomPermeability<dim>::RandomPermeability(std::vector<double> &xi, unsigned int n_terms, double domain_length, double correlation_length, double mu)
  : kl_expansion(n_terms, domain_length, correlation_length, mu)
  , samples(xi)
  {}
  
  template <int dim>
  void RandomField::RandomPermeability<dim>::overwrite_samples(const std::vector<double> &next_sample)
  {
      samples = next_sample;
  }
  
  template <int dim>
  double RandomField::RandomPermeability<dim>::value(const Point<dim>& p)
  {
      return std::exp(kl_expansion.compute_kl_expansion(p, samples));
  }
  
  template class RandomField::RandomPermeability<1>;
  template class RandomField::RandomPermeability<2>;
  

Annotated version of utils/post_processing.py

import re
import matplotlib.pyplot as plt
import numpy as np
import os
 
def plot_from_file(filename):
    if not os.path.exists(filename):
        print(f"Error: {filename} not found. Did you run the C++ code first or is it perhaps in the wrong folder?")
        return
 
    levels, means, variances, samples = [], [], [], []
    
    pattern = re.compile(r"FINAL_STATS Level:(\d+) Mean:([\d\.e+-]+) Var:([\d\.e+-]+) Samples:(\d+)")
 
    with open(filename, 'r') as f:
        for line in f:
            match = pattern.search(line)
            if match:
                levels.append(int(match.group(1)))
                means.append(abs(float(match.group(2)))) # Absolute for log-plot
                variances.append(float(match.group(3)))
                samples.append(int(match.group(4)))
 
    if not levels:
        print("No valid MLMC data found in the file.")
        return
 
    # Convert to arrays for plotting
    levels = np.array(levels)
    
    fig, axs = plt.subplots(1, 3, figsize=(16, 5))
    
    # 1. Variance Decay
    axs[0].plot(levels, variances, 'o-', color='firebrick', label=r'Var[@f$P_l - P_{l-1}@f$]')
    axs[0].set_yscale('log')
    axs[0].set_title('Variance Decay', fontsize=12, fontweight='bold')
    axs[0].set_xlabel('Level')
    axs[0].grid(True, which='both', alpha=0.3)
    axs[0].legend()
 
    # 2. Mean Difference (Bias)
    axs[1].plot(levels, means, 's-', color='royalblue', label=r'|@f$E[P_l - P_{l-1}]@f$|')
    axs[1].set_yscale('log')
    axs[1].set_title('Mean Difference (Bias)', fontsize=12, fontweight='bold')
    axs[1].set_xlabel('Level')
    axs[1].grid(True, which='both', alpha=0.3)
    axs[1].legend()
 
    # 3. Samples per Level
    axs[2].bar(levels, samples, color='seagreen', alpha=0.7)
    axs[2].set_yscale('log')
    axs[2].set_title('Samples per Level (Workload)', fontsize=12, fontweight='bold')
    axs[2].set_xlabel('Level')
    axs[2].set_ylabel('@f$N_l@f$')
    axs[2].grid(axis='y', alpha=0.3, linestyle='--')
 
    plt.tight_layout()
    plt.savefig('mlmc_plots.png', dpi=150) 
    print("Plot saved as 'mlmc_plots.png'")
    plt.show()
 
if __name__ == "__main__":
    plot_from_file('mlmc_results.txt')