add random darcy flow using MLMC

MLMC_random_darcy_flow/CMakeLists.txt

@@ -0,0 +1,32 @@
+# Set the name of the project and target:
+SET(TARGET "mlmc_for_random_darcy_flow")
+
+# Declare all source files the target consists of:
+SET(TARGET_SRC
+  source/KL_expansion.cc
+  source/mlmc.cc
+  source/random_darcy.cc
+  source/random_permeability.cc
+  main.cc
+)
+
+PROJECT(${TARGET} CXX)
+
+
+CMAKE_MINIMUM_REQUIRED(VERSION 3.28.0)
+
+FIND_PACKAGE(deal.II 9.7.1
+  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+  )
+IF(NOT ${deal.II_FOUND})
+  MESSAGE(FATAL_ERROR "\n"
+    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+    )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()
+

MLMC_random_darcy_flow/README.md

@@ -0,0 +1,39 @@
+# 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](https://en.wikipedia.org/wiki/Kardar%E2%80%93Parisi%E2%80%93Zhang_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 @f[ \mathcal{O}(N^{-\frac{1}{2}}) @f], 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 @f[ \mathcal{E}[Q_M]@f], we can use the linearity of expectation to expand the estimator as a telescoping sum: @f[ \mathcal{E}[Q_M] = \mathcal{E}[Q_0] + \sum_{l=1}^{M} \mathcal{E}[Q_l - Q_{l-1}] @f].
+
+One run on level @f[ M-1 @f] is significantly cheaper than one run on level @f[ M @f]. Furthermore, the difference between @f[ Q_M @f] and @f[ Q_{M-1} @f] 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 @f[(\Omega, \mathcal{F}, \mathbb{P})@f] and the domain be @f[D \subset \mathbb{R}^d, d=1,2,3 @f]. We model the hydraulic conductivity as a random field @f[ k = k(x, \omega)@f] on @f[D \times \Omega @f] with a defined mean and covariance. The governing equation is:
+@f[
+-\nabla \cdot ( k(x,\omega)\nabla p(x, \omega)) = f(x).
+@f]
+We assume the right-hand side @f[f(x) @f] is deterministic. Solving this general form is challenging; therefore, we assume @f[ k(x, \omega) @f] follows a log-normal distribution. This involves replacing the conductivity tensor with a scalar-valued field whose logarithm is Gaussian, which guarantees that @f[ k > 0 @f] 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](./doc/images/convergence_mlmc.png)
+
+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](./doc/images/realization.png)
+
+## 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

MLMC_random_darcy_flow/doc/author

@@ -0,0 +1,2 @@
+Jonas Plank <jonas.plank@tum.de>
+

MLMC_random_darcy_flow/doc/builds-on

@@ -0,0 +1,2 @@
+step-5
+

MLMC_random_darcy_flow/doc/entry-name

@@ -0,0 +1,2 @@
+multilevel monte carlo for random darcy flow
+

MLMC_random_darcy_flow/doc/tooltip

@@ -0,0 +1,2 @@
+Implementation of a multilevel monte carlo method for random darcy flow
+

MLMC_random_darcy_flow/include/KL_expansion.h

@@ -0,0 +1,66 @@
+/* -----------------------------------------------------------------------------
+ *
+ * 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_;
+    };
+}

MLMC_random_darcy_flow/include/mlmc.h

@@ -0,0 +1,43 @@
+/* -----------------------------------------------------------------------------
+ *
+ * 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;
+    };
+}
+

MLMC_random_darcy_flow/include/random_darcy.h

@@ -0,0 +1,95 @@
+/* -----------------------------------------------------------------------------
+ *
+ * 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;
+
+    /**
+     * This class is essentially a modification of Step-5. As described in the README,
+     * we treat the finite element solver as a black box. This class demonstrates 
+     * how easily this can be achieved; the only required modifications are 
+     * input/output adjustments and the interface to the random field.
+     */
+    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); 
+
+        /** 
+         * Using an iterative solver for MLMC is generally not ideal.
+         * The matrices can become highly ill-conditioned due to the high 
+         * variations in the random permeability field.
+         */
+        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);                             
+
+        // 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;
+    };
+}

MLMC_random_darcy_flow/include/random_permeability.h

@@ -0,0 +1,37 @@
+/* -----------------------------------------------------------------------------
+ *
+ * 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;
+    };
+}
+