MFEM  v4.0 Finite element discretization library
joule_solver.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
10 // Software Foundation) version 2.1 dated February 1999.
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 static const double SOLVER_TOL = 1.0e-9;
31 static const int SOLVER_MAX_IT = 1000;
32 static int SOLVER_PRINT_LEVEL = 0;
33 static int STATIC_COND = 0;
34
35 // These are defined in joule.cpp
36 void edot_bc(const Vector &x, Vector &E);
37 void e_exact(const Vector &x, double t, Vector &E);
38 void b_exact(const Vector &x, double t, Vector &B);
39 double p_bc(const Vector &x, double t);
40 double t_exact(Vector &x);
41
42 // A Coefficient is an object with a function Eval that returns a double. A
43 // MeshDependentCoefficient returns a different value depending upon the given
44 // mesh attribute, i.e. a "material property".
45 // Somewhat inefficiently, this is achieved using a GridFunction.
47 {
48 private:
49  std::map<int, double> *materialMap;
50  double scaleFactor;
51 public:
52  MeshDependentCoefficient(const std::map<int, double> &inputMap,
53  double scale = 1.0);
55  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
56  void SetScaleFactor(const double &scale) { scaleFactor = scale; }
58  {
59  if (materialMap != NULL) { delete materialMap; }
60  }
61 };
62
63 // This Coefficient is a product of a GridFunction and MeshDependentCoefficient
64 // for example if T (temperature) is a GridFunction and c (heat capacity) is a
65 // MeshDependentCoefficient, this function can compute c*T.
67 {
68 private:
70 public:
72  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
73  void SetMDC(const MeshDependentCoefficient &input_mdc) { mdc = input_mdc; }
74  virtual ~ScaledGFCoefficient() {}
75 };
76
77 /**
78  After spatial discretization, the magnetic diffusion equation can be written
79  as a system of ODEs:
80
81  S0(sigma) P = 0
82  dE/dt = - (M1(sigma) + dt S1(1/mu))^{-1}*(S1(1/mu)*E + sigma Grad P)
83  dB/dt = - Curl(E)
84  dF/dt = (M2(c/k) + dt S2(1))^{-1} (-S2(1) F + Div J)
85  dcT/dt = -Div F + J
86
87  where
88
89  P is the 0-form electrostatic potential,
90  E is the 1-form electric field,
91  B is the 2-form magnetic flux,
92  F is the 2-form thermal flux,
93  T is the 3-form temperature.
94  M is the mass matrix,
95  S is the stiffness matrix,
96  Curl is the curl matrix,
97  Div is the divergence matrix,
98  sigma is the electrical conductivity,
99  k is the thermal conductivity,
100  c is the heat capacity,
101  J is a function of the Joule heating sigma (E dot E).
102
103  Class MagneticDiffusionEOperator represents the right-hand side of the above
104  system of ODEs.
105 */
107 {
108 protected:
109  // These ParFiniteElementSpace objects provide degree-of-freedom mappings.
110  // To create these you must provide the mesh and the definition of the FE
111  // space. These objects are used to create hypre vectors to store the DOFs,
112  // they are used to create grid functions to perform FEM interpolation, and
113  // they are used by bilinear forms.
118
119  // ParBilinearForms are used to create sparse matrices representing discrete
120  // linear operators.
121  ParBilinearForm *a0, *a1, *a2, *m1, *m2, *m3, *s1, *s2;
124
125  // Hypre matrices and vectors for 1-form systems A1 X1 = B1 and 2-form
126  // systems A2 = X2 = B2
127  HypreParMatrix *A0, *A1, *A2, *M1, *M2, *M3;
128  Vector *X0, *X1, *X2, *B0, *B1, *B2, *B3;
129
130  // temporary work vectors
132
133  // HypreSolver is derived from Solver, which is derived from Operator. So a
134  // HypreSolver object has a Mult() operator, which is actually the solver
135  // operation y = A^-1 x i.e. multiplication by A^-1.
136  // HyprePCG is a wrapper for hypre's preconditioned conjugate gradient.
137  mutable HypreSolver * amg_a0;
138  mutable HyprePCG * pcg_a0;
140  mutable HyprePCG * pcg_a2;
141  mutable HypreSolver * ams_a1;
142  mutable HyprePCG * pcg_a1;
143  mutable HypreSolver * dsp_m3;
144  mutable HyprePCG * pcg_m3;
145  mutable HypreSolver * dsp_m1;
146  mutable HyprePCG * pcg_m1;
147  mutable HypreSolver * dsp_m2;
148  mutable HyprePCG * pcg_m2;
149
156
158  double mu, dt_A1, dt_A2;
159
160  // The method builA2 creates the ParBilinearForm a2, the HypreParMatrix A2,
161  // and the solver and preconditioner pcg_a2 and amg_a2. The other build
162  // functions do similar things.
164  void buildA1(double muInv, MeshDependentCoefficient &sigma, double dt);
166  MeshDependentCoefficient &InvTcap, double dt);
169  void buildM3(MeshDependentCoefficient &Tcap);
170  void buildS1(double muInv);
173  void buildCurl(double muInv);
175
176 public:
178  ParFiniteElementSpace &L2FES,
179  ParFiniteElementSpace &HCurlFES,
180  ParFiniteElementSpace &HDivFES,
185  double mu,
186  std::map<int, double> sigmaAttMap,
187  std::map<int, double> TcapacityAttMap,
188  std::map<int, double> InvTcapAttMap,
189  std::map<int, double> InvTcondAttMap
190  );
191
192  // Initialize the fields. This is where restart would go to.
193  void Init(Vector &vx);
194
195  // class TimeDependentOperator is derived from Operator, and class Operator
196  // has the virtual function Mult(x,y) which computes y = A x for some matrix
197  // A, or more generally a nonlinear operator y = A(x).
198  virtual void Mult(const Vector &vx, Vector &dvx_dt) const;
199
200  // Solve the Backward-Euler equation: k = f(x + dt*k, t), for the unknown
201  // slope k. This is the only requirement for high-order SDIRK implicit
202  // integration. This is a virtual function of class TimeDependentOperator.
203  virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
204
205  // Compute B^T M2 B, where M2 is the HDiv mass matrix with permeability
206  // coefficient.
207  // double MagneticEnergy(ParGridFunction &B_gf) const;
208
209  // Compute E^T M1 E, where M1 is the HCurl mass matrix with conductivity
210  // coefficient.
211  double ElectricLosses(ParGridFunction &E_gf) const;
212
213  // E is the input, w is the output which is L2 heating.
214  void GetJouleHeating(ParGridFunction &E_gf, ParGridFunction &w_gf) const;
215
216  void SetTime(const double _t);
217
218  // Write all the hypre matrices and vectors to disk.
219  void Debug(const char *basefilename, double time);
220
221  virtual ~MagneticDiffusionEOperator();
222 };
223
224 // A Coefficient is an object with a function Eval that returns a double. The
225 // JouleHeatingCoefficient object will contain a reference to the electric field
226 // grid function, and the conductivity sigma, and returns sigma E dot E at a
227 // point.
229 {
230 private:
231  ParGridFunction &E_gf;
233 public:
235  ParGridFunction &E_gf_)
236  : E_gf(E_gf_), sigma(sigma_) {}
237  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
239 };
240
241 } // namespace electromagnetics
242
243 } // namespace mfem
244
245 #endif // MFEM_USE_MPI
246
247 #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:27
Base abstract class for time dependent operators.
Definition: operator.hpp:162
void e_exact(const Vector &x, double t, Vector &E)
Definition: joule.cpp:735
double t_exact(Vector &x)
Definition: joule.cpp:749
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:755
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:706
Base class Coefficient that may optionally depend on time.
Definition: coefficient.hpp:31
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:670
void edot_bc(const Vector &x, Vector &E)
Definition: joule.cpp:730
const double alpha
Definition: ex15.cpp:339
Vector data type.
Definition: vector.hpp:48
void b_exact(const Vector &x, double t, Vector &B)
Definition: joule.cpp:742
Class for parallel grid function.
Definition: pgridfunc.hpp:32
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:187
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.