MLMC_random_darcy_flow.h

Go to the documentation of this file.

 
 *   #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:
 *   
 * @endcode
 * 
 * computes our coefficients
 * 
 
 * 
 * computes our eigenvalues based on the previously computed frequencies
 * 
 * @code
 *           void compute_lambda_i(); 
 * @endcode
 * 
 * computes norm values such that our eigenfunctions are normed to one.
 * 
 * @code
 *           void compute_alpha_i();  
 * @endcode
 * 
 * compute the frequencies
 * 
 * @code
 *           void compute_omega_i(); 
 *   
 * @endcode
 * 
 * functions we need for the computation of omega_i
 * 
 * @code
 *           double f_even(double sol);     
 *           double f_odd(double sol);       
 *           double grad_f_even(double sol); 
 *           double grad_f_odd(double sol);  
 *   
 * @endcode
 * 
 * newton since equation is nonlinear
 * 
 * @code
 *           void newton_even(unsigned int index);  
 *           void newton_odd(unsigned int index);   
 *   
 * @endcode
 * 
 * number of KL expansion terms
 * 
 * @code
 *           unsigned int n_terms_;
 * @endcode
 * 
 * domain length
 * 
 * @code
 *           double domain_length_;
 * @endcode
 * 
 * correlation length
 * 
 * @code
 *           double correlation_length_;
 * @endcode
 * 
 * important parameters for the construction of the KL-Expansion in x-direction
 * 
 * @code
 *           std::vector<double> omega_i;
 *           std::vector<double> lambda_i;
 *           std::vector<double> alpha_i;
 *   
 * @endcode
 * 
 * constant mean of our field.
 * 
 * @code
 *           double mu_;
 *       };
 *   }
 * @endcode
 
 
<a name="ann-include/mlmc.h"></a>
<h1>Annotated version of include/mlmc.h</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * 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);                                   
 *   
 * @endcode
 * 
 * generates the required number of samples
 * 
 * @code
 *           std::vector<double> generate_samples();                            
 *   
 * @endcode
 * 
 * store the samples for the computation of mean and variance
 * 
 * @code
 *           void add_sample(double rvalue);                                     
 *   
 *           double compute_mean();
 *   
 *           double compute_variance();
 * @endcode
 * 
 * removes all samples after reached convergence on a level
 * 
 * @code
 *           void clear_samples();
 *   
 *           private:
 *           std::mt19937 rng;                       // Mersenne Twister engine
 *           std::normal_distribution<double> dist;  // Normal distribution
 *           unsigned int num_samples;
 *   
 * @endcode
 * 
 * parameters for error computation
 * 
 * @code
 *           std::vector<double> results;
 *       };
 *   }
 *   
 * @endcode
 
 
<a name="ann-include/random_darcy.h"></a>
<h1>Annotated version of include/random_darcy.h</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * 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();                                                                   
 *   
 * @endcode
 * 
 * Allows switching between coarse and fine triangulations 
 * with virtually no code duplication.
 * 
 * @code
 *           void set_tria(bool fine);                                                       
 *   
 *           void generate_mesh(double domain_length);
 *   
 * @endcode
 * 
 * Implementation follows the logic of deal.II tutorial @ref step_5 "step-5".                                         
 * 
 * @code
 *           void setup_system();                                                            
 *   
 *           void assemble_system(RandomField::RandomPermeability<dim>& random_constant); 
 *   
 *           
 *           void solve();
 *   
 * @endcode
 * 
 * If we are on the first level, we only want to refine the fine triangulation.                                                                           
 * 
 * @code
 *           void refine_grid(bool firstRun);  
 *   
 * @endcode
 * 
 * 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.
 * 
 * @code
 *           void output_results(bool firstRun, unsigned int level);                             
 *   
 * @endcode
 * 
 * This is our Quantity of Interest (QoI).
 * It is defined as: @f$K_{eff} = - \int_{\Gamma_{right}} k \frac{\partial p}{\partial x_1} dx_2@f$
 * 
 * @code
 *           double compute_Keff(RandomField::RandomPermeability<dim>& permeability);
 *   
 *           private:
 * @endcode
 * 
 * define two triangulations
 * 
 * @code
 *           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;
 *       };
 *   }
 * @endcode
 
 
<a name="ann-include/random_permeability.h"></a>
<h1>Annotated version of include/random_permeability.h</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * 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);  
 *   
 * @endcode
 * 
 * overwrite the samples
 * 
 * @code
 *           void overwrite_samples(const std::vector<double> &next_sample);  
 *   
 * @endcode
 * 
 * our evaluator 
 * 
 * @code
 *           double value(const Point<dim>& p);  
 *   
 *           private:
 *           KLExpansion<dim> kl_expansion;
 *           std::vector<double> samples;
 *       };
 *   }
 *   
 * @endcode
 
 
<a name="ann-main.cc"></a>
<h1>Annotated version of main.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   #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;
 *       }
 *   
 * @endcode
 * 
 * Run parameters
 * 
 * @code
 *       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;
 *   
 * @endcode
 * 
 * The number of terms chosen for 1D. For higher dimensions, we use oneD_samples^dim.
 * 
 * @code
 *       unsigned int oneD_samples = 20;
 *   
 * @endcode
 * 
 * The number of levels was not chosen arbitrarily. Convergence is usually so fast 
 * that using more than 4–5 levels is rarely necessary.
 * 
 * @code
 *       unsigned int levels = 5;
 *   
 * @endcode
 * 
 * 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.
 * 
 * @code
 *       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; 
 *   
 * @endcode
 * 
 * In general, these values should be determined by a pilot run. 
 * I chose these specific values to keep the implementation simple.
 * 
 * @code
 *       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());
 *               
 * @endcode
 * 
 * coarse run
 * 
 * @code
 *               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);
 *               }
 *               
 *   
 * @endcode
 * 
 * fine run
 * 
 * @code
 *               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;
 *   }
 *   
 * @endcode
 
 
<a name="ann-source/KL_expansion.cc"></a>
<h1>Annotated version of source/KL_expansion.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * 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++)
 *           {
 * @endcode
 * 
 * even
 * 
 * @code
 *               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));
 *               }
 * @endcode
 * 
 * odd
 * 
 * @code
 *               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++)
 *               {
 * @endcode
 * 
 * both even
 * 
 * @code
 *                   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));
 *                   }
 * @endcode
 * 
 * both odd
 * 
 * @code
 *                   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));
 *                   }
 * @endcode
 * 
 * i odd
 * 
 * @code
 *                   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));
 *                   }
 * @endcode
 * 
 * j odd
 * 
 * @code
 *                   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>;
 *   
 *   
 * @endcode
 
 
<a name="ann-source/mlmc.cc"></a>
<h1>Annotated version of source/mlmc.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * 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>;
 *   
 * @endcode
 
 
<a name="ann-source/random_darcy.cc"></a>
<h1>Annotated version of source/random_darcy.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * 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);
 *   
 * @endcode
 * 
 * 1. Left boundary (indicator 0): u = 1.0
 * 
 * @code
 *       VectorTools::interpolate_boundary_values(dof_handler,
 *                                                0,
 *                                                Functions::ConstantFunction<dim>(1.0),
 *                                                constraints);
 *   
 * @endcode
 * 
 * 2. Right boundary (indicator 1): u = 0.0
 * 
 * @code
 *       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>;
 *   
 * @endcode
 
 
<a name="ann-source/random_permeability.cc"></a>
<h1>Annotated version of source/random_permeability.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * 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>;
 *   
 * @endcode
 
 
<a name="ann-utils/post_processing.py"></a>
<h1>Annotated version of utils/post_processing.py</h1>
@code{.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')
 
    @endcode
 
 
*/