diff --git a/MLMC_random_darcy_flow/CMakeLists.txt b/MLMC_random_darcy_flow/CMakeLists.txt new file mode 100644 index 00000000..86aea094 --- /dev/null +++ b/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() + diff --git a/MLMC_random_darcy_flow/README.md b/MLMC_random_darcy_flow/README.md new file mode 100644 index 00000000..265fb8b1 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/MLMC_random_darcy_flow/doc/author b/MLMC_random_darcy_flow/doc/author new file mode 100644 index 00000000..f3df611d --- /dev/null +++ b/MLMC_random_darcy_flow/doc/author @@ -0,0 +1,2 @@ +Jonas Plank + diff --git a/MLMC_random_darcy_flow/doc/builds-on b/MLMC_random_darcy_flow/doc/builds-on new file mode 100644 index 00000000..389e7378 --- /dev/null +++ b/MLMC_random_darcy_flow/doc/builds-on @@ -0,0 +1,2 @@ +step-5 + diff --git a/MLMC_random_darcy_flow/doc/entry-name b/MLMC_random_darcy_flow/doc/entry-name new file mode 100644 index 00000000..d1b0bee5 --- /dev/null +++ b/MLMC_random_darcy_flow/doc/entry-name @@ -0,0 +1,2 @@ +multilevel monte carlo for random darcy flow + diff --git a/MLMC_random_darcy_flow/doc/images/convergence_mlmc.png b/MLMC_random_darcy_flow/doc/images/convergence_mlmc.png new file mode 100644 index 00000000..1d4a65bd Binary files /dev/null and b/MLMC_random_darcy_flow/doc/images/convergence_mlmc.png differ diff --git a/MLMC_random_darcy_flow/doc/images/realization.png b/MLMC_random_darcy_flow/doc/images/realization.png new file mode 100644 index 00000000..82c0e6e0 Binary files /dev/null and b/MLMC_random_darcy_flow/doc/images/realization.png differ diff --git a/MLMC_random_darcy_flow/doc/tooltip b/MLMC_random_darcy_flow/doc/tooltip new file mode 100644 index 00000000..951d72dd --- /dev/null +++ b/MLMC_random_darcy_flow/doc/tooltip @@ -0,0 +1,2 @@ +Implementation of a multilevel monte carlo method for random darcy flow + diff --git a/MLMC_random_darcy_flow/include/KL_expansion.h b/MLMC_random_darcy_flow/include/KL_expansion.h new file mode 100644 index 00000000..61f877c9 --- /dev/null +++ b/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 +#include + +#include +#include + +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 + class KLExpansion + { + public: + KLExpansion(const unsigned int n_terms, double L, double l, double mu); + + double compute_kl_expansion(const Point& p, std::vector &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 omega_i; + std::vector lambda_i; + std::vector alpha_i; + + // constant mean of our field. + double mu_; + }; +} diff --git a/MLMC_random_darcy_flow/include/mlmc.h b/MLMC_random_darcy_flow/include/mlmc.h new file mode 100644 index 00000000..679e11f9 --- /dev/null +++ b/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 +#include + +namespace MultilevelMonteCarlo +{ + template + class MLMC + { + public: + MLMC(unsigned int oneD_samples); + + // generates the required number of samples + std::vector 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 dist; // Normal distribution + unsigned int num_samples; + + // parameters for error computation + std::vector results; + }; +} + diff --git a/MLMC_random_darcy_flow/include/random_darcy.h b/MLMC_random_darcy_flow/include/random_darcy.h new file mode 100644 index 00000000..5583ff85 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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 + 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& 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& permeability); + + private: + // define two triangulations + Triangulation coarse_tria; + Triangulation fine_tria; + + const FE_Q fe; + DoFHandler dof_handler; + + AffineConstraints constraints; + + SparseMatrix system_matrix; + SparsityPattern sparsity_pattern; + + Vector solution; + Vector system_rhs; + }; +} diff --git a/MLMC_random_darcy_flow/include/random_permeability.h b/MLMC_random_darcy_flow/include/random_permeability.h new file mode 100644 index 00000000..771ecfb1 --- /dev/null +++ b/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 + class RandomPermeability : public Function + { + public: + RandomPermeability(std::vector &first_sample, unsigned int n_terms, double domain_length, double correlation_length, double mu); + + // overwrite the samples + void overwrite_samples(const std::vector &next_sample); + + // our evaluator + double value(const Point& p); + + private: + KLExpansion kl_expansion; + std::vector samples; + }; +} + diff --git a/MLMC_random_darcy_flow/main.cc b/MLMC_random_darcy_flow/main.cc new file mode 100644 index 00000000..020c907e --- /dev/null +++ b/MLMC_random_darcy_flow/main.cc @@ -0,0 +1,127 @@ +#include "include/random_permeability.h" +#include "include/random_darcy.h" +#include "include/mlmc.h" + +#include +#include + +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 runs_per_level{20000, 10000, 5000, 2500, 1000, 250}; + + MultilevelMonteCarlo::MLMC<2> mlmc(oneD_samples); + std::vector 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(int i = 0; i<=levels; i++) + { + std::cout << "Level:" << i << std::endl; + for(unsigned int j = 0; 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; +} + diff --git a/MLMC_random_darcy_flow/source/KL_expansion.cc b/MLMC_random_darcy_flow/source/KL_expansion.cc new file mode 100644 index 00000000..bb148c01 --- /dev/null +++ b/MLMC_random_darcy_flow/source/KL_expansion.cc @@ -0,0 +1,224 @@ +/* ----------------------------------------------------------------------------- + * + * 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 +#include + +template +RandomField::KLExpansion::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 +double RandomField::KLExpansion::compute_kl_expansion(const Point& p, std::vector &samples) +{ + double val = mu_; + if constexpr (dim == 1) + { + for(unsigned int i = 0; i +void RandomField::KLExpansion::compute_omega_i() +{ + for(unsigned int i = 1; i<=n_terms_; i++) + { + if(i%2 == 0) + { + newton_even(i); + } + else + { + newton_odd(i); + } + } +} + +template +void RandomField::KLExpansion::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 +void RandomField::KLExpansion::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 +void RandomField::KLExpansion::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 +void RandomField::KLExpansion::newton_odd(unsigned int index) +{ + double x = M_PI / domain_length_ * (static_cast(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 +double RandomField::KLExpansion::f_odd(double x) +{ + return 1.0 / correlation_length_ - x * std::tan(x * domain_length_ / 2.0); +} + +template +double RandomField::KLExpansion::f_even(double x) +{ + return (1.0 / correlation_length_) * std::tan(x * domain_length_ / 2.0) + x; +} + +template +double RandomField::KLExpansion::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 +double RandomField::KLExpansion::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>; + + diff --git a/MLMC_random_darcy_flow/source/mlmc.cc b/MLMC_random_darcy_flow/source/mlmc.cc new file mode 100644 index 00000000..08396350 --- /dev/null +++ b/MLMC_random_darcy_flow/source/mlmc.cc @@ -0,0 +1,73 @@ +/* ----------------------------------------------------------------------------- + * + * 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 + +template +MultilevelMonteCarlo::MLMC::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 +std::vector MultilevelMonteCarlo::MLMC::generate_samples() +{ + std::vector samples(num_samples); + + for (auto &s : samples) + { + s = dist(rng); + } + return samples; +} + +template +void MultilevelMonteCarlo::MLMC::add_sample(double rvalue) +{ + results.push_back(rvalue); +} + +template +double MultilevelMonteCarlo::MLMC::compute_mean() +{ + double mean = 0.0; + for(int i = 0; i +double MultilevelMonteCarlo::MLMC::compute_variance() +{ + double mean = compute_mean(); + double var = 0.0; + for(int i = 0; i +void MultilevelMonteCarlo::MLMC::clear_samples() +{ + results.clear(); +} + +template class MultilevelMonteCarlo::MLMC<1>; +template class MultilevelMonteCarlo::MLMC<2>; + diff --git a/MLMC_random_darcy_flow/source/random_darcy.cc b/MLMC_random_darcy_flow/source/random_darcy.cc new file mode 100644 index 00000000..b41a5204 --- /dev/null +++ b/MLMC_random_darcy_flow/source/random_darcy.cc @@ -0,0 +1,199 @@ +/* ----------------------------------------------------------------------------- + * + * 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 + +template +Discretization::RandomDarcy::RandomDarcy() +: fe(1) +, dof_handler(coarse_tria) +{} + +template +void Discretization::RandomDarcy::set_tria(bool fine) +{ + if(fine) + { + dof_handler.reinit(fine_tria); + } + else + { + dof_handler.reinit(coarse_tria); + } +} + +template +void Discretization::RandomDarcy::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 +void Discretization::RandomDarcy::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); + + // 1. Left boundary (indicator 0): u = 1.0 + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Functions::ConstantFunction(1.0), + constraints); + + // 2. Right boundary (indicator 1): u = 0.0 + VectorTools::interpolate_boundary_values(dof_handler, + 1, + Functions::ConstantFunction(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 +void Discretization::RandomDarcy::assemble_system(RandomField::RandomPermeability& permeability) +{ + const QGauss quadrature_formula(fe.degree + 1); + + FEValues 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 cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + std::vector 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 +void Discretization::RandomDarcy::solve() +{ + SparseDirectUMFPACK A_direct; + + solution = system_rhs; + A_direct.solve(system_matrix, solution); + + constraints.distribute(solution); +} + +template +void Discretization::RandomDarcy::refine_grid(bool firstRun) +{ + fine_tria.refine_global(1); + if(!firstRun) + { + coarse_tria.refine_global(1); + } +} + +template +void Discretization::RandomDarcy::output_results(bool firstRun, unsigned int level) +{ + DataOut 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 +double Discretization::RandomDarcy::compute_Keff(RandomField::RandomPermeability& permeability) +{ + const QGauss face_quadrature_formula(fe.degree + 1); + FEFaceValues fe_face_values(fe, + face_quadrature_formula, + update_gradients | update_normal_vectors | + update_JxW_values | update_quadrature_points); + + std::vector> 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>; + diff --git a/MLMC_random_darcy_flow/source/random_permeability.cc b/MLMC_random_darcy_flow/source/random_permeability.cc new file mode 100644 index 00000000..71317805 --- /dev/null +++ b/MLMC_random_darcy_flow/source/random_permeability.cc @@ -0,0 +1,32 @@ +/* ----------------------------------------------------------------------------- + * + * 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 +RandomField::RandomPermeability::RandomPermeability(std::vector &xi, unsigned int n_terms, double domain_length, double correlation_length, double mu) +: samples(xi) +, kl_expansion(n_terms, domain_length, correlation_length, mu) +{} + +template +void RandomField::RandomPermeability::overwrite_samples(const std::vector &next_sample) +{ + samples = next_sample; +} + +template +double RandomField::RandomPermeability::value(const Point& p) +{ + return std::exp(kl_expansion.compute_kl_expansion(p, samples)); +} + +template class RandomField::RandomPermeability<1>; +template class RandomField::RandomPermeability<2>; + diff --git a/MLMC_random_darcy_flow/utils/post_processing.py b/MLMC_random_darcy_flow/utils/post_processing.py new file mode 100644 index 00000000..4af35f70 --- /dev/null +++ b/MLMC_random_darcy_flow/utils/post_processing.py @@ -0,0 +1,65 @@ +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[$P_l - P_{l-1}$]') + 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'|$E[P_l - P_{l-1}]$|') + 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('$N_l$') + 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') + + \ No newline at end of file