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