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