MFEM  v4.6.0 Finite element discretization library
maxwell.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
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 // ------------------------------------------------------------------
13 // Maxwell Miniapp: Simple Full-Wave Electromagnetic Simulation Code
14 // ------------------------------------------------------------------
15 //
16 // This miniapp solves a simple 3D full-wave electromagnetic problem using the
17 // coupled, first-order equations:
18 //
19 // epsilon dE/dt = Curl 1/mu B - sigma E - J
20 // dB/dt = - Curl E
21 //
22 // The permittivity function is that of the vacuum with an optional dielectric
23 // sphere. The permeability function is that of the vacuum with an optional
24 // diamagnetic or paramagnetic spherical shell. The optional conductivity
25 // function is also a user-defined sphere.
26 //
27 // The optional current density is a pulse of current in the shape of a cylinder
28 // with a time dependence resembling the derivative of a Gaussian distribution.
29 //
30 // Boundary conditions can be 'natural' meaning zero tangential current,
31 // 'Dirichlet' which sets the time-derivative of the tangential components of E,
32 // or 'absorbing' (we use a simple Sommerfeld first order absorbing boundary
33 // condition).
34 //
35 // We discretize the electric field with H(Curl) finite elements (Nedelec edge
36 // elements) and the magnetic flux with H(Div) finite elements (Raviart-Thomas
37 // elements).
38 //
39 // The symplectic time integration algorithm used below is designed to conserve
40 // energy unless lossy materials or absorbing boundary conditions are used.
41 // When losses are expected, the algorithm uses an implicit method which
42 // includes the loss operators in the left hand side of the linear system.
43 //
44 // For increased accuracy the time integration order can be set to 2, 3, or 4
45 // (the default is 1st order).
46 //
47 // Compile with: make maxwell
48 //
49 // Sample runs:
50 //
51 // Current source in a sphere with absorbing boundary conditions:
52 // mpirun -np 4 maxwell -m ../../data/ball-nurbs.mesh -rs 2
53 // -abcs '-1'
54 // -dp '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5'
55 //
56 // Current source in a metal sphere with dielectric and conducting materials:
57 // mpirun -np 4 maxwell -m ../../data/ball-nurbs.mesh -rs 2
58 // -dbcs '-1'
59 // -dp '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5'
60 // -cs '0.0 0.0 -0.5 .2 3e6'
61 // -ds '0.0 0.0 0.5 .2 10'
62 //
63 // Current source in a metal box:
64 // mpirun -np 4 maxwell -m ../../data/fichera.mesh -rs 3
65 // -ts 0.25 -tf 10 -dbcs '-1'
66 // -dp '-0.5 -0.5 0.0 -0.5 -0.5 1.0 0.1 1 .5 1'
67 //
68 // Current source with a mixture of absorbing and reflecting boundaries:
69 // mpirun -np 4 maxwell -m ../../data/fichera.mesh -rs 3
70 // -ts 0.25 -tf 10
71 // -dp '-0.5 -0.5 0.0 -0.5 -0.5 1.0 0.1 1 .5 1'
72 // -dbcs '4 8 19 21' -abcs '5 18'
73 //
74 // By default the sources and fields are all zero:
75 // mpirun -np 4 maxwell
76
77 #include "maxwell_solver.hpp"
78 #include <fstream>
79 #include <iostream>
80
81 using namespace std;
82 using namespace mfem;
83 using namespace mfem::common;
84 using namespace mfem::electromagnetics;
85
86 // Permittivity Function
87 static Vector ds_params_(0); // Center, Radius, and Permittivity
88 // of dielectric sphere
89 double dielectric_sphere(const Vector &);
90 double epsilon(const Vector &x) { return dielectric_sphere(x); }
91
92 // Permeability Function
93 static Vector ms_params_(0); // Center, Inner and Outer Radii, and
94 // Permeability of magnetic shell
95 double magnetic_shell(const Vector &);
96 double muInv(const Vector & x) { return 1.0/magnetic_shell(x); }
97
98 // Conductivity Function
99 static Vector cs_params_(0); // Center, Radius, and Conductivity
100 // of conductive sphere
101 double conductive_sphere(const Vector &);
102 double sigma(const Vector &x) { return conductive_sphere(x); }
103
104 // Current Density Function
105 static Vector dp_params_(0); // Axis Start, Axis End, Rod Radius,
106 // Total Current of Rod, and Frequency
107 void dipole_pulse(const Vector &x, double t, Vector &j);
108 void j_src(const Vector &x, double t, Vector &j) { dipole_pulse(x, t, j); }
109
110 // dE/dt Boundary Condition: The following function returns zero but any time
111 // dependent function could be used.
112 void dEdtBCFunc(const Vector &x, double t, Vector &E);
113
114 // The following functions return zero but they could be modified to set initial
115 // conditions for the electric and magnetic fields
116 void EFieldFunc(const Vector &, Vector&);
117 void BFieldFunc(const Vector &, Vector&);
118
119 // Scale factor between input time units and seconds
120 static double tScale_ = 1e-9; // Input time in nanosecond
121
122 int SnapTimeStep(double tmax, double dtmax, double & dt);
123
124 // Prints the program's logo to the given output stream
125 void display_banner(ostream & os);
126
127 int main(int argc, char *argv[])
128 {
129  Mpi::Init();
130  Hypre::Init();
131
132  if ( Mpi::Root() ) { display_banner(cout); }
133
134  // Parse command-line options.
135  const char *mesh_file = "../../data/ball-nurbs.mesh";
136  int sOrder = 1;
137  int tOrder = 1;
138  int serial_ref_levels = 0;
139  int parallel_ref_levels = 0;
140  bool visualization = true;
141  bool visit = true;
142  double dt = 1.0e-12;
143  double dtsf = 0.95;
144  double ti = 0.0;
145  double ts = 1.0;
146  double tf = 40.0;
147
148  Array<int> abcs;
149  Array<int> dbcs;
150
151  OptionsParser args(argc, argv);
153  "Mesh file to use.");
155  "Finite element order (polynomial degree).");
157  "Time integration order.");
159  "Number of serial refinement levels.");
161  "Number of parallel refinement levels.");
163  "Used to reduce the time step below the upper bound.");
165  "Beginning of time interval to simulate (ns).");
167  "End of time interval to simulate (ns).");
169  "Time between snapshots (ns).");
171  "Center, Radius, and Permittivity of Dielectric Sphere");
174  "of Magnetic Shell");
176  "Center, Radius, and Conductivity of Conductive Sphere");
178  "Axis End Points, Radius, Amplitude, "
179  "Pulse Center (ns), Pulse Width (ns)");
181  "Absorbing Boundary Condition Surfaces");
183  "Dirichlet Boundary Condition Surfaces");
185  "--no-visualization",
186  "Enable or disable GLVis visualization.");
188  "--no-visualization",
189  "Enable or disable VisIt visualization.");
190  args.Parse();
191  if (!args.Good())
192  {
193  if (Mpi::Root())
194  {
195  args.PrintUsage(cout);
196  }
197  return 1;
198  }
199  if (Mpi::Root())
200  {
201  args.PrintOptions(cout);
202  }
203
204  // Read the (serial) mesh from the given mesh file on all processors. We can
205  // handle triangular, quadrilateral, tetrahedral, hexahedral, surface and
206  // volume meshes with the same code.
207  Mesh *mesh;
208  ifstream imesh(mesh_file);
209  if (!imesh)
210  {
211  if (Mpi::Root())
212  {
213  cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
214  }
215  return 2;
216  }
217  mesh = new Mesh(imesh, 1, 1);
218  imesh.close();
219
220  // Project a NURBS mesh to a piecewise-quadratic curved mesh
221  if (mesh->NURBSext)
222  {
223  mesh->UniformRefinement();
224  if (serial_ref_levels > 0) { serial_ref_levels--; }
225
226  mesh->SetCurvature(2);
227  }
228
229  // Refine the serial mesh on all processors to increase the resolution. In
230  // this example we do 'ref_levels' of uniform refinement.
231  for (int l = 0; l < serial_ref_levels; l++)
232  {
233  mesh->UniformRefinement();
234  }
235
236  // Define a parallel mesh by a partitioning of the serial mesh. Refine this
237  // mesh further in parallel to increase the resolution. Once the parallel
238  // mesh is defined, the serial mesh can be deleted.
239  ParMesh pmesh(MPI_COMM_WORLD, *mesh);
240  delete mesh;
241
242  // Refine this mesh in parallel to increase the resolution.
243  for (int l = 0; l < parallel_ref_levels; l++)
244  {
245  pmesh.UniformRefinement();
246  }
247
248  // Create the Electromagnetic solver
249  MaxwellSolver Maxwell(pmesh, sOrder,
250  ( ds_params_.Size() > 0 ) ? epsilon : NULL,
251  ( ms_params_.Size() > 0 ) ? muInv : NULL,
252  ( cs_params_.Size() > 0 ) ? sigma : NULL,
253  ( dp_params_.Size() > 0 ) ? j_src : NULL,
254  abcs, dbcs,
255  ( dbcs.Size() > 0 ) ? dEdtBCFunc : NULL
256  );
257
258  // Display the current number of DoFs in each finite element space
259  Maxwell.PrintSizes();
260
261  // Set the initial conditions for both the electric and magnetic fields
264
265  Maxwell.SetInitialEField(EFieldCoef);
266  Maxwell.SetInitialBField(BFieldCoef);
267
268  // Compute the energy of the initial fields
269  double energy = Maxwell.GetEnergy();
270  if ( Mpi::Root() )
271  {
272  cout << "Energy(" << ti << "ns): " << energy << "J" << endl;
273  }
274
275  // Approximate the largest stable time step
276  double dtmax = Maxwell.GetMaximumTimeStep();
277
278  // Convert times from nanoseconds to seconds
279  ti *= tScale_;
280  tf *= tScale_;
281  ts *= tScale_;
282
283  if ( Mpi::Root() )
284  {
285  cout << "Maximum Time Step: " << dtmax / tScale_ << "ns" << endl;
286  }
287
288  // Round down the time step so that tf-ti is an integer multiple of dt
289  int nsteps = SnapTimeStep(tf-ti, dtsf * dtmax, dt);
290  if ( Mpi::Root() )
291  {
292  cout << "Number of Time Steps: " << nsteps << endl;
293  cout << "Time Step Size: " << dt / tScale_ << "ns" << endl;
294  }
295
296  // Create the ODE solver
297  SIAVSolver siaSolver(tOrder);
298  siaSolver.Init(Maxwell.GetNegCurl(), Maxwell);
299
300
301  // Initialize GLVis visualization
302  if (visualization)
303  {
304  Maxwell.InitializeGLVis();
305  }
306
307  // Initialize VisIt visualization
308  VisItDataCollection visit_dc("Maxwell-Parallel", &pmesh);
309
310  double t = ti;
311  Maxwell.SetTime(t);
312
313  if ( visit )
314  {
315  Maxwell.RegisterVisItFields(visit_dc);
316  }
317
318  // Write initial fields to disk for VisIt
319  if ( visit )
320  {
321  Maxwell.WriteVisItFields(0);
322  }
323
324  // Send the initial condition by socket to a GLVis server.
325  if (visualization)
326  {
327  Maxwell.DisplayToGLVis();
328  }
329
330  // The main time evolution loop.
331  int it = 1;
332  while (t < tf)
333  {
334  // Run the simulation until a snapshot is needed
335  siaSolver.Run(Maxwell.GetBField(), Maxwell.GetEField(), t, dt,
336  max(t + dt, ti + ts * it));
337
338  // Approximate the current energy if the fields
339  energy = Maxwell.GetEnergy();
340  if ( Mpi::Root() )
341  {
342  cout << "Energy(" << t/tScale_ << "ns): " << energy << "J" << endl;
343  }
344
345  // Update local DoFs with current true DoFs
346  Maxwell.SyncGridFuncs();
347
348  // Write fields to disk for VisIt
349  if ( visit )
350  {
351  Maxwell.WriteVisItFields(it);
352  }
353
354  // Send the solution by socket to a GLVis server.
355  if (visualization)
356  {
357  Maxwell.DisplayToGLVis();
358  }
359
360  it++;
361  }
362
363  return 0;
364 }
365
366 // Print the Maxwell ascii logo to the given ostream
367 void display_banner(ostream & os)
368 {
369  os << " ___ ____ " << endl
370  << " / | / / __ __ " << endl
371  << " / |_/ _ /__ ___ _____ _ __ ____ | | | | " << endl
372  << " / \\__ \\ \\ \\/ /\\ \\/ \\/ // __ \\| | | | "
373  << endl
374  << " / /|_/ // __ \\_> < \\ /\\ ___/| |_| |__ " << endl
375  << "/___/ /_ /(____ /__/\\_ \\ \\/\\_/ \\___ >____/____/ " << endl
376  << " \\/ \\/ \\/ \\/ " << endl
377  << flush;
378 }
379
380 // A sphere with constant permittivity. The sphere has a radius, center, and
381 // permittivity specified on the command line and stored in ds_params_.
382 double dielectric_sphere(const Vector &x)
383 {
384  double r2 = 0.0;
385
386  for (int i=0; i<x.Size(); i++)
387  {
388  r2 += (x(i)-ds_params_(i))*(x(i)-ds_params_(i));
389  }
390
391  if ( sqrt(r2) <= ds_params_(x.Size()) )
392  {
393  return ds_params_(x.Size()+1) * epsilon0_;
394  }
395  return epsilon0_;
396 }
397
398 // A spherical shell with constant permeability. The sphere has inner and outer
399 // radii, center, and relative permeability specified on the command line and
400 // stored in ms_params_.
401 double magnetic_shell(const Vector &x)
402 {
403  double r2 = 0.0;
404
405  for (int i=0; i<x.Size(); i++)
406  {
407  r2 += (x(i)-ms_params_(i))*(x(i)-ms_params_(i));
408  }
409
410  if ( sqrt(r2) >= ms_params_(x.Size()) &&
411  sqrt(r2) <= ms_params_(x.Size()+1) )
412  {
413  return mu0_*ms_params_(x.Size()+2);
414  }
415  return mu0_;
416 }
417
418 // A sphere with constant conductivity. The sphere has a radius, center, and
419 // conductivity specified on the command line and stored in ls_params_.
420 double conductive_sphere(const Vector &x)
421 {
422  double r2 = 0.0;
423
424  for (int i=0; i<x.Size(); i++)
425  {
426  r2 += (x(i)-cs_params_(i))*(x(i)-cs_params_(i));
427  }
428
429  if ( sqrt(r2) <= cs_params_(x.Size()) )
430  {
431  return cs_params_(x.Size()+1);
432  }
433  return 0.0;
434 }
435
436 // A cylindrical rod of current density. The rod has two axis end points, a
437 // radius, a current amplitude in Amperes, a center time, and a width. All of
438 // these parameters are stored in dp_params_.
439 void dipole_pulse(const Vector &x, double t, Vector &j)
440 {
441  MFEM_ASSERT(x.Size() == 3, "current source requires 3D space.");
442
443  j.SetSize(x.Size());
444  j = 0.0;
445
446  Vector v(x.Size()); // Normalized Axis vector
447  Vector xu(x.Size()); // x vector relative to the axis end-point
448
449  xu = x;
450
451  for (int i=0; i<x.Size(); i++)
452  {
453  xu[i] -= dp_params_[i];
454  v[i] = dp_params_[x.Size()+i] - dp_params_[i];
455  }
456
457  double h = v.Norml2();
458
459  if ( h == 0.0 )
460  {
461  return;
462  }
463  v /= h;
464
465  double r = dp_params_[2*x.Size()+0];
466  double a = dp_params_[2*x.Size()+1] * tScale_;
467  double b = dp_params_[2*x.Size()+2] * tScale_;
468  double c = dp_params_[2*x.Size()+3] * tScale_;
469
470  double xv = xu * v;
471
472  // Compute perpendicular vector from axis to x
474
475  double xp = xu.Norml2();
476
477  if ( xv >= 0.0 && xv <= h && xp <= r )
478  {
479  j = v;
480  }
481
482  j *= a * (t - b) * exp(-0.5 * pow((t-b)/c, 2.0)) / (c * c);
483 }
484
485 void
486 EFieldFunc(const Vector &x, Vector &E)
487 {
488  E.SetSize(3);
489  E = 0.0;
490 }
491
492 void
493 BFieldFunc(const Vector &x, Vector &B)
494 {
495  B.SetSize(3);
496  B = 0.0;
497 }
498
499 void
500 dEdtBCFunc(const Vector &x, double t, Vector &dE)
501 {
502  dE.SetSize(3);
503  dE = 0.0;
504 }
505
506 int
507 SnapTimeStep(double tmax, double dtmax, double & dt)
508 {
509  double dsteps = tmax/dtmax;
510
511  int nsteps = static_cast<int>(pow(10,(int)ceil(log10(dsteps))));
512
513  for (int i=1; i<=5; i++)
514  {
515  int a = (int)ceil(log10(dsteps/pow(5.0,i)));
516  int nstepsi = (int)pow(5,i)*max(1,(int)pow(10,a));
517
518  nsteps = min(nsteps,nstepsi);
519  }
520
521  dt = tmax / nsteps;
522
523  return nsteps;
524 }
virtual void Init(Operator &P, TimeDependentOperator &F)
Definition: ode.cpp:916
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
int SnapTimeStep(double tmax, double dtmax, double &dt)
Definition: maxwell.cpp:507
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
virtual void SetTime(const double t_)
Set the current time.
Definition: operator.hpp:360
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
double muInv(const Vector &x)
Definition: maxwell.cpp:96
STL namespace.
void SetInitialBField(VectorCoefficient &BFieldCoef)
double magnetic_shell(const Vector &)
Definition: maxwell.cpp:401
void BFieldFunc(const Vector &, Vector &)
Definition: maxwell.cpp:493
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
void RegisterVisItFields(VisItDataCollection &visit_dc)
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
Data collection with VisIt I/O routines.
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition: mesh.cpp:5635
int main(int argc, char *argv[])
Definition: maxwell.cpp:127
void display_banner(ostream &os)
Definition: maxwell.cpp:367
A general vector function coefficient.
double dielectric_sphere(const Vector &)
Definition: maxwell.cpp:382
void dipole_pulse(const Vector &x, double t, Vector &j)
Definition: maxwell.cpp:439
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
void dEdtBCFunc(const Vector &x, double t, Vector &E)
Definition: maxwell.cpp:500
void EFieldFunc(const Vector &, Vector &)
Definition: maxwell.cpp:486
double conductive_sphere(const Vector &)
Definition: maxwell.cpp:420
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:277
double epsilon
Definition: maxwell.cpp:126
void j_src(const Vector &x, double t, Vector &j)
Definition: maxwell.cpp:108
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:835
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
RefCoord t[3]
virtual void Run(Vector &q, Vector &p, double &t, double &dt, double tf)
Definition: ode.hpp:579
Vector data type.
Definition: vector.hpp:58
double sigma(const Vector &x)
Definition: maxwell.cpp:102
Class for parallel meshes.
Definition: pmesh.hpp:32
void SetInitialEField(VectorCoefficient &EFieldCoef)
Variable order Symplectic Integration Algorithm (orders 1-4)
Definition: ode.hpp:611