MFEM  v3.4
Finite element discretization library
maxwell.cpp
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
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
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::miniapps;
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 // depenent 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_Session mpi(argc, argv);
130 
131  if ( mpi.Root() ) { display_banner(cout); }
132 
133  // Parse command-line options.
134  const char *mesh_file = "../../data/ball-nurbs.mesh";
135  int sOrder = 1;
136  int tOrder = 1;
137  int serial_ref_levels = 0;
138  int parallel_ref_levels = 0;
139  bool visualization = true;
140  bool visit = true;
141  double dt = 1.0e-12;
142  double dtsf = 0.95;
143  double ti = 0.0;
144  double ts = 1.0;
145  double tf = 40.0;
146 
147  Array<int> abcs;
148  Array<int> dbcs;
149 
150  OptionsParser args(argc, argv);
151  args.AddOption(&mesh_file, "-m", "--mesh",
152  "Mesh file to use.");
153  args.AddOption(&sOrder, "-so", "--spatial-order",
154  "Finite element order (polynomial degree).");
155  args.AddOption(&tOrder, "-to", "--temporal-order",
156  "Time integration order.");
157  args.AddOption(&serial_ref_levels, "-rs", "--serial-ref-levels",
158  "Number of serial refinement levels.");
159  args.AddOption(&parallel_ref_levels, "-rp", "--parallel-ref-levels",
160  "Number of parallel refinement levels.");
161  args.AddOption(&dtsf, "-sf", "--dt-safety-factor",
162  "Used to reduce the time step below the upper bound.");
163  args.AddOption(&ti, "-ti", "--initial-time",
164  "Beginning of time interval to simulate (ns).");
165  args.AddOption(&tf, "-tf", "--final-time",
166  "End of time interval to simulate (ns).");
167  args.AddOption(&ts, "-ts", "--snapshot-time",
168  "Time between snapshots (ns).");
169  args.AddOption(&ds_params_, "-ds", "--dielectric-sphere-params",
170  "Center, Radius, and Permittivity of Dielectric Sphere");
171  args.AddOption(&ms_params_, "-ms", "--magnetic-shell-params",
172  "Center, Inner Radius, Outer Radius, and Permeability "
173  "of Magnetic Shell");
174  args.AddOption(&cs_params_, "-cs", "--conductive-sphere-params",
175  "Center, Radius, and Conductivity of Conductive Sphere");
176  args.AddOption(&dp_params_, "-dp", "--dipole-pulse-params",
177  "Axis End Points, Radius, Amplitude, "
178  "Pulse Center (ns), Pulse Width (ns)");
179  args.AddOption(&abcs, "-abcs", "--absorbing-bc-surf",
180  "Absorbing Boundary Condition Surfaces");
181  args.AddOption(&dbcs, "-dbcs", "--dirichlet-bc-surf",
182  "Dirichlet Boundary Condition Surfaces");
183  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
184  "--no-visualization",
185  "Enable or disable GLVis visualization.");
186  args.AddOption(&visit, "-visit", "--visit", "-no-visit",
187  "--no-visualization",
188  "Enable or disable VisIt visualization.");
189  args.Parse();
190  if (!args.Good())
191  {
192  if (mpi.Root())
193  {
194  args.PrintUsage(cout);
195  }
196  return 1;
197  }
198  if (mpi.Root())
199  {
200  args.PrintOptions(cout);
201  }
202 
203  // Read the (serial) mesh from the given mesh file on all processors. We can
204  // handle triangular, quadrilateral, tetrahedral, hexahedral, surface and
205  // volume meshes with the same code.
206  Mesh *mesh;
207  ifstream imesh(mesh_file);
208  if (!imesh)
209  {
210  if (mpi.Root())
211  {
212  cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
213  }
214  return 2;
215  }
216  mesh = new Mesh(imesh, 1, 1);
217  imesh.close();
218 
219  // Project a NURBS mesh to a piecewise-quadratic curved mesh
220  if (mesh->NURBSext)
221  {
222  mesh->UniformRefinement();
223  if (serial_ref_levels > 0) { serial_ref_levels--; }
224 
225  mesh->SetCurvature(2);
226  }
227 
228  // Refine the serial mesh on all processors to increase the resolution. In
229  // this example we do 'ref_levels' of uniform refinement.
230  for (int l = 0; l < serial_ref_levels; l++)
231  {
232  mesh->UniformRefinement();
233  }
234 
235  // Define a parallel mesh by a partitioning of the serial mesh. Refine this
236  // mesh further in parallel to increase the resolution. Once the parallel
237  // mesh is defined, the serial mesh can be deleted.
238  ParMesh pmesh(MPI_COMM_WORLD, *mesh);
239  delete mesh;
240 
241  // Refine this mesh in parallel to increase the resolution.
242  for (int l = 0; l < parallel_ref_levels; l++)
243  {
244  pmesh.UniformRefinement();
245  }
246 
247  // Create the Electromagnetic solver
248  MaxwellSolver Maxwell(pmesh, sOrder,
249  ( ds_params_.Size() > 0 ) ? epsilon : NULL,
250  ( ms_params_.Size() > 0 ) ? muInv : NULL,
251  ( cs_params_.Size() > 0 ) ? sigma : NULL,
252  ( dp_params_.Size() > 0 ) ? j_src : NULL,
253  abcs, dbcs,
254  ( dbcs.Size() > 0 ) ? dEdtBCFunc : NULL
255  );
256 
257  // Display the current number of DoFs in each finite element space
258  Maxwell.PrintSizes();
259 
260  // Set the initial conditions for both the electric and magnetic fields
263 
264  Maxwell.SetInitialEField(EFieldCoef);
265  Maxwell.SetInitialBField(BFieldCoef);
266 
267  // Compute the energy of the initial fields
268  double energy = Maxwell.GetEnergy();
269  if ( mpi.Root() )
270  {
271  cout << "Energy(" << ti << "ns): " << energy << "J" << endl;
272  }
273 
274  // Approximate the largest stable time step
275  double dtmax = Maxwell.GetMaximumTimeStep();
276 
277  // Convert times from nanoseconds to seconds
278  ti *= tScale_;
279  tf *= tScale_;
280  ts *= tScale_;
281 
282  if ( mpi.Root() )
283  {
284  cout << "Maximum Time Step: " << dtmax / tScale_ << "ns" << endl;
285  }
286 
287  // Round down the time step so that tf-ti is an integer multiple of dt
288  int nsteps = SnapTimeStep(tf-ti, dtsf * dtmax, dt);
289  if ( mpi.Root() )
290  {
291  cout << "Number of Time Steps: " << nsteps << endl;
292  cout << "Time Step Size: " << dt / tScale_ << "ns" << endl;
293  }
294 
295  // Create the ODE solver
296  SIAVSolver siaSolver(tOrder);
297  siaSolver.Init(Maxwell.GetNegCurl(), Maxwell);
298 
299 
300  // Initialize GLVis visualization
301  if (visualization)
302  {
303  Maxwell.InitializeGLVis();
304  }
305 
306  // Initialize VisIt visualization
307  VisItDataCollection visit_dc("Maxwell-Parallel", &pmesh);
308 
309  double t = ti;
310  Maxwell.SetTime(t);
311 
312  if ( visit )
313  {
314  Maxwell.RegisterVisItFields(visit_dc);
315  }
316 
317  // Write initial fields to disk for VisIt
318  if ( visit )
319  {
320  Maxwell.WriteVisItFields(0);
321  }
322 
323  // Send the initial condition by socket to a GLVis server.
324  if (visualization)
325  {
326  Maxwell.DisplayToGLVis();
327  }
328 
329  // The main time evolution loop.
330  int it = 1;
331  while (t < tf)
332  {
333  // Run the simulation until a snapshot is needed
334  siaSolver.Run(Maxwell.GetBField(), Maxwell.GetEField(), t, dt,
335  max(t + dt, ti + ts * it));
336 
337  // Approximate the current energy if the fields
338  energy = Maxwell.GetEnergy();
339  if ( mpi.Root() )
340  {
341  cout << "Energy(" << t/tScale_ << "ns): " << energy << "J" << endl;
342  }
343 
344  // Update local DoFs with current true DoFs
345  Maxwell.SyncGridFuncs();
346 
347  // Write fields to disk for VisIt
348  if ( visit )
349  {
350  Maxwell.WriteVisItFields(it);
351  }
352 
353  // Send the solution by socket to a GLVis server.
354  if (visualization)
355  {
356  Maxwell.DisplayToGLVis();
357  }
358 
359  it++;
360  }
361 
362  return 0;
363 }
364 
365 // Print the Maxwell ascii logo to the given ostream
366 void display_banner(ostream & os)
367 {
368  os << " ___ ____ " << endl
369  << " / | / / __ __ " << endl
370  << " / |_/ _ /__ ___ _____ _ __ ____ | | | | " << endl
371  << " / \\__ \\ \\ \\/ /\\ \\/ \\/ // __ \\| | | | "
372  << endl
373  << " / /|_/ // __ \\_> < \\ /\\ ___/| |_| |__ " << endl
374  << "/___/ /_ /(____ /__/\\_ \\ \\/\\_/ \\___ >____/____/ " << endl
375  << " \\/ \\/ \\/ \\/ " << endl
376  << flush;
377 }
378 
379 // A sphere with constant permittivity. The sphere has a radius, center, and
380 // permittivity specified on the command line and stored in ds_params_.
381 double dielectric_sphere(const Vector &x)
382 {
383  double r2 = 0.0;
384 
385  for (int i=0; i<x.Size(); i++)
386  {
387  r2 += (x(i)-ds_params_(i))*(x(i)-ds_params_(i));
388  }
389 
390  if ( sqrt(r2) <= ds_params_(x.Size()) )
391  {
392  return ds_params_(x.Size()+1) * epsilon0_;
393  }
394  return epsilon0_;
395 }
396 
397 // A spherical shell with constant permeability. The sphere has inner and outer
398 // radii, center, and relative permeability specified on the command line and
399 // stored in ms_params_.
400 double magnetic_shell(const Vector &x)
401 {
402  double r2 = 0.0;
403 
404  for (int i=0; i<x.Size(); i++)
405  {
406  r2 += (x(i)-ms_params_(i))*(x(i)-ms_params_(i));
407  }
408 
409  if ( sqrt(r2) >= ms_params_(x.Size()) &&
410  sqrt(r2) <= ms_params_(x.Size()+1) )
411  {
412  return mu0_*ms_params_(x.Size()+2);
413  }
414  return mu0_;
415 }
416 
417 // A sphere with constant conductivity. The sphere has a radius, center, and
418 // conductivity specified on the command line and stored in ls_params_.
419 double conductive_sphere(const Vector &x)
420 {
421  double r2 = 0.0;
422 
423  for (int i=0; i<x.Size(); i++)
424  {
425  r2 += (x(i)-cs_params_(i))*(x(i)-cs_params_(i));
426  }
427 
428  if ( sqrt(r2) <= cs_params_(x.Size()) )
429  {
430  return cs_params_(x.Size()+1);
431  }
432  return 0.0;
433 }
434 
435 // A cylindrical rod of current density. The rod has two axis end points, a
436 // radus, a current amplitude in Amperes, a center time, and a width. All of
437 // these parameters are stored in dp_params_.
438 void dipole_pulse(const Vector &x, double t, Vector &j)
439 {
440  MFEM_ASSERT(x.Size() == 3, "current source requires 3D space.");
441 
442  j.SetSize(x.Size());
443  j = 0.0;
444 
445  Vector v(x.Size()); // Normalized Axis vector
446  Vector xu(x.Size()); // x vector relative to the axis end-point
447 
448  xu = x;
449 
450  for (int i=0; i<x.Size(); i++)
451  {
452  xu[i] -= dp_params_[i];
453  v[i] = dp_params_[x.Size()+i] - dp_params_[i];
454  }
455 
456  double h = v.Norml2();
457 
458  if ( h == 0.0 )
459  {
460  return;
461  }
462  v /= h;
463 
464  double r = dp_params_[2*x.Size()+0];
465  double a = dp_params_[2*x.Size()+1] * tScale_;
466  double b = dp_params_[2*x.Size()+2] * tScale_;
467  double c = dp_params_[2*x.Size()+3] * tScale_;
468 
469  double xv = xu * v;
470 
471  // Compute perpendicular vector from axis to x
472  xu.Add(-xv, v);
473 
474  double xp = xu.Norml2();
475 
476  if ( xv >= 0.0 && xv <= h && xp <= r )
477  {
478  j = v;
479  }
480 
481  j *= a * (t - b) * exp(-0.5 * pow((t-b)/c, 2)) / (c * c);
482 }
483 
484 void
485 EFieldFunc(const Vector &x, Vector &E)
486 {
487  E.SetSize(3);
488  E = 0.0;
489 }
490 
491 void
492 BFieldFunc(const Vector &x, Vector &B)
493 {
494  B.SetSize(3);
495  B = 0.0;
496 }
497 
498 void
499 dEdtBCFunc(const Vector &x, double t, Vector &dE)
500 {
501  dE.SetSize(3);
502  dE = 0.0;
503 }
504 
505 int
506 SnapTimeStep(double tmax, double dtmax, double & dt)
507 {
508  double dsteps = tmax/dtmax;
509 
510  int nsteps = pow(10,(int)ceil(log10(dsteps)));
511 
512  for (int i=1; i<=5; i++)
513  {
514  int a = (int)ceil(log10(dsteps/pow(5.0,i)));
515  int nstepsi = (int)pow(5,i)*max(1,(int)pow(10,a));
516 
517  nsteps = min(nsteps,nstepsi);
518  }
519 
520  dt = tmax / nsteps;
521 
522  return nsteps;
523 }
virtual void Init(Operator &P, TimeDependentOperator &F)
Definition: ode.cpp:565
bool Root() const
Return true if WorldRank() == 0.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:328
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
double dielectric_sphere(const Vector &)
Definition: maxwell.cpp:381
int Size() const
Returns the size of the vector.
Definition: vector.hpp:120
virtual void SetTime(const double _t)
Set the current time.
Definition: operator.hpp:181
void EFieldFunc(const Vector &, Vector &)
Definition: maxwell.cpp:485
bool Good() const
Definition: optparser.hpp:120
void display_banner(ostream &os)
Definition: maxwell.cpp:366
STL namespace.
void SetInitialBField(VectorCoefficient &BFieldCoef)
A simple convenience class that calls MPI_Init() at construction and MPI_Finalize() at destruction...
void RegisterVisItFields(VisItDataCollection &visit_dc)
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6741
void dEdtBCFunc(const Vector &x, double t, Vector &E)
Definition: maxwell.cpp:499
Data collection with VisIt I/O routines.
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3355
double muInv(const Vector &x)
Definition: maxwell.cpp:96
void dipole_pulse(const Vector &x, double t, Vector &j)
Definition: maxwell.cpp:438
int SnapTimeStep(double tmax, double dtmax, double &dt)
Definition: maxwell.cpp:506
int main(int argc, char *argv[])
Definition: maxwell.cpp:127
void j_src(const Vector &x, double t, Vector &j)
Definition: maxwell.cpp:108
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)
Definition: optparser.hpp:74
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:176
double conductive_sphere(const Vector &)
Definition: maxwell.cpp:419
double magnetic_shell(const Vector &)
Definition: maxwell.cpp:400
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:666
int Size() const
Logical size of the array.
Definition: array.hpp:133
double epsilon(const Vector &x)
Definition: maxwell.cpp:90
virtual void Run(Vector &q, Vector &p, double &t, double &dt, double tf)
Definition: ode.hpp:332
Vector data type.
Definition: vector.hpp:48
void BFieldFunc(const Vector &, Vector &)
Definition: maxwell.cpp:492
Class for parallel meshes.
Definition: pmesh.hpp:32
void SetInitialEField(VectorCoefficient &EFieldCoef)
double sigma(const Vector &x)
Definition: maxwell.cpp:102