MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
spde_solver.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details
11
12#ifndef SPDE_SOLVERS_HPP
13#define SPDE_SOLVERS_HPP
14
15#include <ostream>
16#include <unordered_map>
17#include "mfem.hpp"
18
19namespace mfem
20{
21namespace spde
22{
23
25
27{
28 Boundary() = default;
29
30 /// Print the information specifying the boundary conditions.
31 void PrintInfo(std::ostream &os = mfem::out) const;
32
33 /// Verify that all defined boundaries are actually defined on the mesh, i.e.
34 /// if the keys of boundary attributes appear in the boundary attributes of
35 /// the mesh.
36 void VerifyDefinedBoundaries(const Mesh &mesh) const;
37
38 /// Computes the error for each defined boundary attribute by calling back to
39 /// the IntegrateBC.
40 void ComputeBoundaryError(const ParGridFunction &solution);
41
42 /// Helper function to compute the coefficients alpha, beta, gamma (see in
43 /// `IntegrateBC`) for a given boundary attribute.
45 real_t &gamma);
46
47 /// Add a homogeneous boundary condition to the boundary.
48 void AddHomogeneousBoundaryCondition(int boundary, BoundaryType type);
49
50 /// Add a inhomogeneous Dirichlet boundary condition to the boundary.
52 real_t coefficient);
53
54 /// Set the robin coefficient for the boundary.
55 void SetRobinCoefficient(real_t coefficient);
56
57 /// Map to assign homogeneous boundary conditions to defined boundary types.
58 std::map<int, BoundaryType> boundary_attributes;
59 /// Coefficient for inhomogeneous Dirichlet boundary conditions.
60 std::map<int, real_t> dirichlet_coefficients;
61 /// Coefficient for Robin boundary conditions (n.grad(u) + coeff u = 0) on
62 /// defined boundaries.
64};
65
66/// IntegrateBC function from ex27p.cpp. For boundary verification.
67/// Compute the average value of alpha*n.Grad(sol) + beta*sol over the boundary
68/// attributes marked in bdr_marker. Also computes the L2 norm of
69/// alpha*n.Grad(sol) + beta*sol - gamma over the same boundary.
70real_t IntegrateBC(const ParGridFunction &x, const Array<int> &bdr,
71 real_t alpha, real_t beta, real_t gamma, real_t &glb_err);
72
73/// Solver for the SPDE method based on a rational approximation with the AAA
74/// algorithm. The SPDE method is described in the paper
75/// Lindgren, F., Rue, H., Lindström, J. (2011). An explicit link between
76/// Gaussian fields and Gaussian Markov random fields: the stochastic partial
77/// differential equation approach. Journal of the Royal Statistical Society:
78/// Series B (Statistical Methodology), 73(4), 423–498.
79/// https://doi.org/10.1111/j.1467-9868.2011.00777.x
80///
81/// The solver solves the SPDE problem defined as
82/// (A)^-alpha u = b
83/// where A is
84/// A = div ( Theta(x) grad + Id ) u(x)
85/// and alpha is given as
86/// alpha = (2 nu + dim) / 2.
87/// Theta (anisotropy tensor) and nu (smoothness) can be specified in the
88/// constructor. Traditionally, the SPDE method requires the specification of
89/// a white noise right hands side. SPDESolver accepts arbitrary right hand
90/// sides but the solver has only been tested with white noise.
92{
93public:
94 /// Constructor.
95 /// @param nu The coefficient nu, smoothness of the solution.
96 /// @param bc Boundary conditions.
97 /// @param fespace Finite element space.
98 /// @param l1 Correlation length in x
99 /// @param l2 Correlation length in y
100 /// @param l3 Correlation length in z
101 /// @param e1 Rotation angle in x
102 /// @param e2 Rotation angle in y
103 /// @param e3 Rotation angle in z
104 SPDESolver(real_t nu, const Boundary &bc, ParFiniteElementSpace *fespace,
105 real_t l1 = 0.1, real_t l2 = 0.1, real_t l3 = 0.1,
106 real_t e1 = 0.0, real_t e2 = 0.0, real_t e3 = 0.0);
107
108 /// Destructor.
109 ~SPDESolver();
110
111 /// Solve the SPDE for a given right hand side b. May alter b if
112 /// the exponent (alpha) is larger than 1. We avoid copying by default. If
113 /// you need b later on, make a copy of it before calling this function.
115
116 /// Set up the random field generator. If called more than once it resets the
117 /// generator.
118 void SetupRandomFieldGenerator(int seed=0);
119
120 /// Generate a random field. Calls back to solve but generates the stochastic
121 /// load b internally.
123
124 /// Construct the normalization coefficient eta of the white noise right hands
125 /// side.
127 real_t l2, real_t l3,
128 int dim);
129
130 /// Construct the second order tensor (matrix coefficient) Theta from the
131 /// equation R^T(e1,e2,e3) diag(l1,l2,l3) R (e1,e2,e3).
133 real_t e1, real_t e2, real_t e3,
134 real_t nu, int dim);
135
136 /// Set the print level
137 void SetPrintLevel(int print_level) {print_level_ = print_level;}
138
139private:
140 /// The rational approximation of the SPDE results in multiple
141 /// reaction-diffusion PDEs that need to be solved. This call solves the PDE
142 /// (div Theta grad + alpha I)^exponent x = beta b.
144 real_t beta, int exponent = 1);
145
146 /// Lift the solution to satisfy the inhomogeneous boundary conditions.
147 void LiftSolution(ParGridFunction &x);
148
149 // Each PDE gives rise to a linear system. This call solves the linear system
150 // with PCG and Boomer AMG preconditioner.
151 void SolveLinearSystem(const HypreParMatrix *Op);
152
153 /// Activate repeated solve capabilities. E.g. if the PDE is of the form
154 /// A^N x = b. This method solves the PDE A x = b for the first time, and
155 /// then uses the solution as RHS for the next solve and so forth.
156 void ActivateRepeatedSolve() { repeated_solve_ = true; }
157
158 /// Single solve only.
159 void DeactivateRepeatedSolve() { repeated_solve_ = false; }
160
161 /// Writes the solution of the PDE from the previous call to Solve() to the
162 /// linear from b (with appropriate transformations).
163 void UpdateRHS(ParLinearForm &b) const;
164
165 // Compute the coefficients for the rational approximation of the solution.
166 void ComputeRationalCoefficients(real_t exponent);
167
168 // Bilinear forms and corresponding matrices for the solver.
169 ParBilinearForm k_;
170 ParBilinearForm m_;
171 HypreParMatrix stiffness_;
172 HypreParMatrix mass_bc_;
173
174 // Transformation matrices (needed to construct the linear systems and
175 // solutions)
176 const SparseMatrix *restriction_matrix_ = nullptr;
177 const Operator *prolongation_matrix_ = nullptr;
178
179 // Members to solve the linear system.
180 Vector X_;
181 Vector B_;
182
183 // Information of the finite element space.
184 Array<int> ess_tdof_list_;
185 ParFiniteElementSpace *fespace_ptr_;
186
187 // Boundary conditions
188 const Boundary &bc_;
189 Array<int> dbc_marker_; // Markers for Dirichlet boundary conditions.
190 Array<int> rbc_marker_; // Markers for Robin boundary conditions.
191
192 // Coefficients for the rational approximation of the solution.
193 Array<real_t> coeffs_;
194 Array<real_t> poles_;
195
196 // Exponents of the operator
197 real_t nu_ = 0.0;
198 real_t alpha_ = 0.0;
199 int integer_order_of_exponent_ = 0;
200
201 // Correlation length
202 real_t l1_ = 0.1;
203 real_t l2_ = 0.1;
204 real_t l3_ = 0.1;
205 real_t e1_ = 0.0;
206 real_t e2_ = 0.0;
207 real_t e3_ = 0.0;
208
209 // Print level
210 int print_level_ = 1;
211
212 // Member to switch to repeated solve capabilities.
213 bool repeated_solve_ = false;
214 bool integer_order_ = false;
215 bool apply_lift_ = false;
216
217 WhiteGaussianNoiseDomainLFIntegrator *integ = nullptr;
218 ParLinearForm * b_wn = nullptr;
219};
220
221} // namespace spde
222} // namespace mfem
223
224#endif // SPDE_SOLVERS_HPP
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
Mesh data type.
Definition mesh.hpp:56
Abstract parallel finite element space.
Definition pfespace.hpp:29
Class for parallel grid function.
Definition pgridfunc.hpp:33
Class for parallel linear form.
void Solve(ParLinearForm &b, ParGridFunction &x)
static DenseMatrix ConstructMatrixCoefficient(real_t l1, real_t l2, real_t l3, real_t e1, real_t e2, real_t e3, real_t nu, int dim)
static real_t ConstructNormalizationCoefficient(real_t nu, real_t l1, real_t l2, real_t l3, int dim)
void GenerateRandomField(ParGridFunction &x)
void SetupRandomFieldGenerator(int seed=0)
void SetPrintLevel(int print_level)
Set the print level.
SPDESolver(real_t nu, const Boundary &bc, ParFiniteElementSpace *fespace, real_t l1=0.1, real_t l2=0.1, real_t l3=0.1, real_t e1=0.0, real_t e2=0.0, real_t e3=0.0)
Vector beta
const real_t alpha
Definition ex15.cpp:369
int dim
Definition ex24.cpp:53
real_t b
Definition lissajous.cpp:42
real_t IntegrateBC(const ParGridFunction &x, const Array< int > &bdr, real_t alpha, real_t beta, real_t gamma, real_t &glb_err)
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition globals.hpp:66
float real_t
Definition config.hpp:43
std::map< int, BoundaryType > boundary_attributes
Map to assign homogeneous boundary conditions to defined boundary types.
std::map< int, real_t > dirichlet_coefficients
Coefficient for inhomogeneous Dirichlet boundary conditions.
void PrintInfo(std::ostream &os=mfem::out) const
Print the information specifying the boundary conditions.
void ComputeBoundaryError(const ParGridFunction &solution)
void AddInhomogeneousDirichletBoundaryCondition(int boundary, real_t coefficient)
Add a inhomogeneous Dirichlet boundary condition to the boundary.
void SetRobinCoefficient(real_t coefficient)
Set the robin coefficient for the boundary.
void UpdateIntegrationCoefficients(int i, real_t &alpha, real_t &beta, real_t &gamma)
void VerifyDefinedBoundaries(const Mesh &mesh) const
void AddHomogeneousBoundaryCondition(int boundary, BoundaryType type)
Add a homogeneous boundary condition to the boundary.