MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
joule_solver.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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 MFEM_JOULE_SOLVER
13 #define MFEM_JOULE_SOLVER
14 
15 #include "../common/pfem_extras.hpp"
16 
17 #ifdef MFEM_USE_MPI
18 
19 #include <memory>
20 #include <iostream>
21 #include <fstream>
22 
23 namespace mfem
24 {
25 
26 namespace electromagnetics
27 {
28 
29 // Some global variable for convenience
30 const double SOLVER_TOL = 1.0e-9;
31 const int SOLVER_MAX_IT = 1000;
32 // Initialized in joule.cpp and used in joule_solver.cpp:
33 extern int SOLVER_PRINT_LEVEL;
34 extern int STATIC_COND;
35 
36 // These are defined in joule.cpp
37 void edot_bc(const Vector &x, Vector &E);
38 void e_exact(const Vector &x, double t, Vector &E);
39 void b_exact(const Vector &x, double t, Vector &B);
40 double p_bc(const Vector &x, double t);
41 double t_exact(const Vector &x);
42 
43 // A Coefficient is an object with a function Eval that returns a double. A
44 // MeshDependentCoefficient returns a different value depending upon the given
45 // mesh attribute, i.e. a "material property".
46 // Somewhat inefficiently, this is achieved using a GridFunction.
48 {
49 private:
50  std::map<int, double> *materialMap;
51  double scaleFactor;
52 public:
53  MeshDependentCoefficient(const std::map<int, double> &inputMap,
54  double scale = 1.0);
56  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
57  void SetScaleFactor(const double &scale) { scaleFactor = scale; }
59  {
60  if (materialMap != NULL) { delete materialMap; }
61  }
62 };
63 
64 // This Coefficient is a product of a GridFunction and MeshDependentCoefficient
65 // for example if T (temperature) is a GridFunction and c (heat capacity) is a
66 // MeshDependentCoefficient, this function can compute c*T.
68 {
69 private:
71 public:
73  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
74  void SetMDC(const MeshDependentCoefficient &input_mdc) { mdc = input_mdc; }
75  virtual ~ScaledGFCoefficient() {}
76 };
77 
78 /**
79  After spatial discretization, the magnetic diffusion equation can be written
80  as a system of ODEs:
81 
82  S0(sigma) P = 0
83  dE/dt = - (M1(sigma) + dt S1(1/mu))^{-1}*(S1(1/mu)*E + sigma Grad P)
84  dB/dt = - Curl(E)
85  dF/dt = (M2(c/k) + dt S2(1))^{-1} (-S2(1) F + Div J)
86  dcT/dt = -Div F + J
87 
88  where
89 
90  P is the 0-form electrostatic potential,
91  E is the 1-form electric field,
92  B is the 2-form magnetic flux,
93  F is the 2-form thermal flux,
94  T is the 3-form temperature.
95  M is the mass matrix,
96  S is the stiffness matrix,
97  Curl is the curl matrix,
98  Div is the divergence matrix,
99  sigma is the electrical conductivity,
100  k is the thermal conductivity,
101  c is the heat capacity,
102  J is a function of the Joule heating sigma (E dot E).
103 
104  Class MagneticDiffusionEOperator represents the right-hand side of the above
105  system of ODEs.
106 */
108 {
109 protected:
110  // These ParFiniteElementSpace objects provide degree-of-freedom mappings.
111  // To create these you must provide the mesh and the definition of the FE
112  // space. These objects are used to create hypre vectors to store the DOFs,
113  // they are used to create grid functions to perform FEM interpolation, and
114  // they are used by bilinear forms.
119 
120  // ParBilinearForms are used to create sparse matrices representing discrete
121  // linear operators.
122  ParBilinearForm *a0, *a1, *a2, *m1, *m2, *m3, *s1, *s2;
125 
126  // Hypre matrices and vectors for 1-form systems A1 X1 = B1 and 2-form
127  // systems A2 = X2 = B2
128  HypreParMatrix *A0, *A1, *A2, *M1, *M2, *M3;
129  Vector *X0, *X1, *X2, *B0, *B1, *B2, *B3;
130 
131  // temporary work vectors
133 
134  // HypreSolver is derived from Solver, which is derived from Operator. So a
135  // HypreSolver object has a Mult() operator, which is actually the solver
136  // operation y = A^-1 x i.e. multiplication by A^-1.
137  // HyprePCG is a wrapper for hypre's preconditioned conjugate gradient.
138  mutable HypreSolver * amg_a0;
139  mutable HyprePCG * pcg_a0;
140  mutable HypreSolver * ads_a2;
141  mutable HyprePCG * pcg_a2;
142  mutable HypreSolver * ams_a1;
143  mutable HyprePCG * pcg_a1;
144  mutable HypreSolver * dsp_m3;
145  mutable HyprePCG * pcg_m3;
146  mutable HypreSolver * dsp_m1;
147  mutable HyprePCG * pcg_m1;
148  mutable HypreSolver * dsp_m2;
149  mutable HyprePCG * pcg_m2;
150 
157 
159  double mu, dt_A1, dt_A2;
160 
161  // The method builA2 creates the ParBilinearForm a2, the HypreParMatrix A2,
162  // and the solver and preconditioner pcg_a2 and amg_a2. The other build
163  // functions do similar things.
165  void buildA1(double muInv, MeshDependentCoefficient &sigma, double dt);
167  MeshDependentCoefficient &InvTcap, double dt);
170  void buildM3(MeshDependentCoefficient &Tcap);
171  void buildS1(double muInv);
173  void buildGrad();
174  void buildCurl(double muInv);
176 
177 public:
179  ParFiniteElementSpace &L2FES,
180  ParFiniteElementSpace &HCurlFES,
181  ParFiniteElementSpace &HDivFES,
182  ParFiniteElementSpace &HGradFES,
186  double mu,
187  std::map<int, double> sigmaAttMap,
188  std::map<int, double> TcapacityAttMap,
189  std::map<int, double> InvTcapAttMap,
190  std::map<int, double> InvTcondAttMap
191  );
192 
193  // Initialize the fields. This is where restart would go to.
194  void Init(Vector &vx);
195 
196  // class TimeDependentOperator is derived from Operator, and class Operator
197  // has the virtual function Mult(x,y) which computes y = A x for some matrix
198  // A, or more generally a nonlinear operator y = A(x).
199  virtual void Mult(const Vector &vx, Vector &dvx_dt) const;
200 
201  // Solve the Backward-Euler equation: k = f(x + dt*k, t), for the unknown
202  // slope k. This is the only requirement for high-order SDIRK implicit
203  // integration. This is a virtual function of class TimeDependentOperator.
204  virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
205 
206  // Compute B^T M2 B, where M2 is the HDiv mass matrix with permeability
207  // coefficient.
208  // double MagneticEnergy(ParGridFunction &B_gf) const;
209 
210  // Compute E^T M1 E, where M1 is the HCurl mass matrix with conductivity
211  // coefficient.
212  double ElectricLosses(ParGridFunction &E_gf) const;
213 
214  // E is the input, w is the output which is L2 heating.
215  void GetJouleHeating(ParGridFunction &E_gf, ParGridFunction &w_gf) const;
216 
217  void SetTime(const double _t);
218 
219  // Write all the hypre matrices and vectors to disk.
220  void Debug(const char *basefilename, double time);
221 
222  virtual ~MagneticDiffusionEOperator();
223 };
224 
225 // A Coefficient is an object with a function Eval that returns a double. The
226 // JouleHeatingCoefficient object will contain a reference to the electric field
227 // grid function, and the conductivity sigma, and returns sigma E dot E at a
228 // point.
230 {
231 private:
232  ParGridFunction &E_gf;
234 public:
236  ParGridFunction &E_gf_)
237  : E_gf(E_gf_), sigma(sigma_) {}
238  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
240 };
241 
242 } // namespace electromagnetics
243 
244 } // namespace mfem
245 
246 #endif // MFEM_USE_MPI
247 
248 #endif // MFEM_JOULE_SOLVER
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
void buildDiv(MeshDependentCoefficient &InvTcap)
void buildM2(MeshDependentCoefficient &alpha)
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
Base abstract class for first order time dependent operators.
Definition: operator.hpp:267
void e_exact(const Vector &x, double t, Vector &E)
Definition: joule.cpp:740
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k)
Solve the equation: k = f(x + dt k, t), for the unknown k at the current time t.
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
void Debug(const char *basefilename, double time)
void GetJouleHeating(ParGridFunction &E_gf, ParGridFunction &w_gf) const
void buildA1(double muInv, MeshDependentCoefficient &sigma, double dt)
ScaledGFCoefficient(GridFunction *gf, MeshDependentCoefficient &input_mdc)
double p_bc(const Vector &x, double t)
Definition: joule.cpp:760
void buildS2(MeshDependentCoefficient &alpha)
void buildM3(MeshDependentCoefficient &Tcap)
void buildM1(MeshDependentCoefficient &sigma)
MagneticDiffusionEOperator(int len, ParFiniteElementSpace &L2FES, ParFiniteElementSpace &HCurlFES, ParFiniteElementSpace &HDivFES, ParFiniteElementSpace &HGradFES, Array< int > &ess_bdr, Array< int > &thermal_ess_bdr, Array< int > &poisson_ess_bdr, double mu, std::map< int, double > sigmaAttMap, std::map< int, double > TcapacityAttMap, std::map< int, double > InvTcapAttMap, std::map< int, double > InvTcondAttMap)
double muInv(const Vector &x)
Definition: maxwell.cpp:96
double ElectricLosses(ParGridFunction &E_gf) const
JouleHeatingCoefficient(const MeshDependentCoefficient &sigma_, ParGridFunction &E_gf_)
MeshDependentCoefficient(const std::map< int, double > &inputMap, double scale=1.0)
PCG solver in hypre.
Definition: hypre.hpp:759
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:39
void buildA2(MeshDependentCoefficient &InvTcond, MeshDependentCoefficient &InvTcap, double dt)
Class for integration point with weight.
Definition: intrules.hpp:25
Class for parallel bilinear form using different test and trial FE spaces.
void SetMDC(const MeshDependentCoefficient &input_mdc)
Class for parallel bilinear form.
Abstract class for hypre&#39;s solvers and preconditioners.
Definition: hypre.hpp:700
void edot_bc(const Vector &x, Vector &E)
Definition: joule.cpp:735
const double alpha
Definition: ex15.cpp:336
Vector data type.
Definition: vector.hpp:51
double t_exact(const Vector &x)
Definition: joule.cpp:754
void b_exact(const Vector &x, double t, Vector &B)
Definition: joule.cpp:747
Class for parallel grid function.
Definition: pgridfunc.hpp:32
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:181
void buildA0(MeshDependentCoefficient &sigma)
virtual void Mult(const Vector &vx, Vector &dvx_dt) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
void SetTime(const double _t)
Set the current time.