MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
maxwell_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_MAXWELL_SOLVER
13 #define MFEM_MAXWELL_SOLVER
14 
15 #include "../common/pfem_extras.hpp"
16 #include "../common/mesh_extras.hpp"
17 #include "electromagnetics.hpp"
18 
19 #ifdef MFEM_USE_MPI
20 
21 #include <string>
22 #include <map>
23 
24 using namespace std;
25 using namespace mfem;
26 
27 namespace mfem
28 {
29 
33 
34 namespace electromagnetics
35 {
36 
38 {
39 public:
40  MaxwellSolver(ParMesh & pmesh, int sOrder,
41  double (*eps )(const Vector&),
42  double (*muInv )(const Vector&),
43  double (*sigma )(const Vector&),
44  void (*j_src )(const Vector&, double, Vector&),
45  Array<int> & abcs, Array<int> & dbcs,
46  void (*dEdt_bc )(const Vector&, double, Vector&));
47 
48  ~MaxwellSolver();
49 
50  int GetLogging() const { return logging_; }
51  void SetLogging(int logging) { logging_ = logging; }
52 
53  HYPRE_Int GetProblemSize();
54 
55  void PrintSizes();
56 
57  void SetInitialEField(VectorCoefficient & EFieldCoef);
58  void SetInitialBField(VectorCoefficient & BFieldCoef);
59 
60  void Mult(const Vector &B, Vector &dEdt) const;
61 
62  void ImplicitSolve(const double dt, const Vector &x, Vector &k);
63 
64  double GetMaximumTimeStep() const;
65 
66  double GetEnergy() const;
67 
68  Operator & GetNegCurl() { return *NegCurl_; }
69 
70  Vector & GetEField() { return *E_; }
71  Vector & GetBField() { return *B_; }
72 
73  void SyncGridFuncs();
74 
75  void RegisterVisItFields(VisItDataCollection & visit_dc);
76 
77  void WriteVisItFields(int it = 0);
78 
79  void InitializeGLVis();
80 
81  void DisplayToGLVis();
82 
83 private:
84 
85  // This method alters mutable member data
86  void setupSolver(const int idt, const double dt) const;
87 
88  void implicitSolve(const double dt, const Vector &x, Vector &k) const;
89 
90  int myid_;
91  int num_procs_;
92  int order_;
93  int logging_;
94 
95  bool lossy_;
96 
97  double dtMax_; // Maximum stable time step
98  double dtScale_; // Used to scale dt before converting to an integer
99 
100  ParMesh * pmesh_;
101 
102  ND_ParFESpace * HCurlFESpace_;
103  RT_ParFESpace * HDivFESpace_;
104 
105  ParBilinearForm * hDivMassMuInv_;
106  ParBilinearForm * hCurlLosses_;
107  ParMixedBilinearForm * weakCurlMuInv_;
108 
109  ParDiscreteCurlOperator * Curl_;
110 
111  ParGridFunction * e_; // Electric Field (HCurl)
112  ParGridFunction * b_; // Magnetic Flux (HDiv)
113  ParGridFunction * j_; // Volumetric Current Density (HCurl)
114  ParGridFunction * dedt_; // Time Derivative of Electric Field (HCurl)
115  ParGridFunction * rhs_; // Dual of displacement current, rhs vector (HCurl)
116  ParLinearForm * jd_; // Dual of current density (HCurl)
117 
118  HypreParMatrix * M1Losses_;
119  HypreParMatrix * M2MuInv_;
120  HypreParMatrix * NegCurl_;
121  HypreParMatrix * WeakCurlMuInv_;
122  HypreParVector * E_; // Current value of the electric field DoFs
123  HypreParVector * B_; // Current value of the magnetic flux DoFs
124  mutable HypreParVector * HD_; // Used in energy calculation
125  mutable HypreParVector * RHS_;
126 
127  Coefficient * epsCoef_; // Electric Permittivity Coefficient
128  Coefficient * muInvCoef_; // Magnetic Permeability Coefficient
129  Coefficient * sigmaCoef_; // Electric Conductivity Coefficient
130  Coefficient * etaInvCoef_; // Admittance Coefficient
131  VectorCoefficient * eCoef_; // Initial Electric Field
132  VectorCoefficient * bCoef_; // Initial Magnetic Flux
133  VectorCoefficient * jCoef_; // Time dependent current density
134  VectorCoefficient * dEdtBCCoef_; // Time dependent boundary condition
135 
136  double (*eps_ )(const Vector&);
137  double (*muInv_ )(const Vector&);
138  double (*sigma_ )(const Vector&);
139  void (*j_src_ )(const Vector&, double, Vector&);
140 
141  // Array of 0's and 1's marking the location of absorbing surfaces
142  Array<int> abc_marker_;
143 
144  // Array of 0's and 1's marking the location of Dirichlet boundaries
145  Array<int> dbc_marker_;
146  void (*dEdt_bc_)(const Vector&, double, Vector&);
147 
148  // Dirichlet degrees of freedom
149  Array<int> dbc_dofs_;
150 
151  // High order symplectic integration requires partial time steps of differing
152  // lengths. If losses are present the system matrix includes a portion scaled
153  // by the time step. Consequently, high order time integration requires
154  // different system matrices. The following maps contain various objects that
155  // depend on the time step.
156  mutable std::map<int, ParBilinearForm *> a1_;
157  mutable std::map<int, HypreParMatrix *> A1_;
158  mutable std::map<int, Coefficient *> dtCoef_;
159  mutable std::map<int, Coefficient *> dtSigmaCoef_;
160  mutable std::map<int, Coefficient *> dtEtaInvCoef_;
161  mutable std::map<int, HypreDiagScale *> diagScale_;
162  mutable std::map<int, HyprePCG *> pcg_;
163 
164  // Data collection used to write VisIt files
165  VisItDataCollection * visit_dc_;
166 
167  // Sockets used to communicate with GLVis
168  std::map<std::string, socketstream*> socks_;
169 };
170 
171 } // namespace electromagnetics
172 
173 } // namespace mfem
174 
175 #endif // MFEM_USE_MPI
176 
177 #endif // MFEM_MAXWELL_SOLVER
Base class for vector Coefficients that optionally depend on time and space.
Base abstract class for first order time dependent operators.
Definition: operator.hpp:267
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
Class for parallel linear form.
Definition: plinearform.hpp:26
Data collection with VisIt I/O routines.
double muInv(const Vector &x)
Definition: maxwell.cpp:96
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:70
void j_src(const Vector &x, double t, Vector &j)
Definition: maxwell.cpp:108
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:39
Class for parallel bilinear form using different test and trial FE spaces.
Class for parallel bilinear form.
Vector data type.
Definition: vector.hpp:51
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:181
Class for parallel meshes.
Definition: pmesh.hpp:32
double sigma(const Vector &x)
Definition: maxwell.cpp:102