MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
volta.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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// Volta Miniapp: Simple Electrostatics Simulation Code
14// -----------------------------------------------------
15//
16// This miniapp solves a simple 2D or 3D electrostatic problem.
17//
18// Div eps Grad Phi = rho
19//
20// The permittivity function is that of the vacuum with an optional dielectric
21// sphere. The charge density is either zero or a user defined sphere of charge.
22//
23// Boundary conditions for the electric potential consist of a user defined
24// piecewise constant potential or a potential leading to a user selected
25// uniform electric field.
26//
27// We discretize the electric potential with H1 finite elements. The electric
28// field E is discretized with Nedelec finite elements.
29//
30// Compile with: make volta
31//
32// Sample runs:
33//
34// Three point charges within a large metal enclosure:
35// mpirun -np 4 volta -m ../../data/inline-quad.mesh
36// -pc '0.5 0.42 20 0.5 0.5 -12 0.5 0.545 15'
37// -dbcs '1 2 3 4' -dbcv '0 0 0 0'
38//
39// A cylinder at constant voltage in a square, grounded metal pipe:
40// mpirun -np 4 volta -m ../../data/square-disc.mesh
41// -dbcs '1 2 3 4 5 6 7 8' -dbcv '0 0 0 0 1 1 1 1'
42//
43// A cylinder with a constant surface charge density in a square,
44// grounded metal pipe:
45// mpirun -np 4 volta -m ../../data/square-disc.mesh
46// -nbcs '5 6 7 8' -nbcv '5e-11 5e-11 5e-11 5e-11'
47// -dbcs '1 2 3 4'
48//
49// A cylindrical voltaic pile within a grounded metal sphere:
50// mpirun -np 4 volta -dbcs 1 -vp '0 -0.5 0 0 0.5 0 0.2 1'
51//
52// A charged sphere, off-center, within a grounded metal sphere:
53// mpirun -np 4 volta -dbcs 1 -cs '0.0 0.5 0.0 0.2 2.0e-11'
54//
55// A dielectric sphere suspended in a uniform electric field:
56// mpirun -np 4 volta -dbcs 1 -dbcg -ds '0.0 0.0 0.0 0.2 8.0'
57//
58// An example using piecewise constant permittivity values
59// mpirun -np 4 volta -m llnl.mesh -dbcs '4' -dbcv '0'
60// -cs '8.5 8.5 17 1.57' -pwe '1 1 1 0.001'
61//
62// By default the sources and fields are all zero:
63// mpirun -np 4 volta
64
65#include "volta_solver.hpp"
66#include <fstream>
67#include <iostream>
68
69using namespace std;
70using namespace mfem;
71using namespace mfem::electromagnetics;
72
73// Permittivity Functions
75
76static Vector pw_eps_(0); // Piecewise permittivity values
77static Vector ds_params_(0); // Center, Radius, and Permittivity
78// of dielectric sphere
80
81// Charge Density Function
82static Vector cs_params_(0); // Center, Radius, and Total Charge
83// of charged sphere
85
86// Point Charges
87static Vector pc_params_(0); // Point charge locations and charges
88
89// Polarization
90static Vector vp_params_(0); // Axis Start, Axis End, Cylinder Radius,
91// and Polarization Magnitude
92void voltaic_pile(const Vector &, Vector &);
93
94// Phi Boundary Condition
95static Vector e_uniform_(0);
97
98// Prints the program's logo to the given output stream
99void display_banner(ostream & os);
100
101int main(int argc, char *argv[])
102{
103 Mpi::Init(argc, argv);
104 Hypre::Init();
105
106 if ( Mpi::Root() ) { display_banner(cout); }
107
108 // Parse command-line options.
109 const char *mesh_file = "../../data/ball-nurbs.mesh";
110 int order = 1;
111 int maxit = 100;
112 int serial_ref_levels = 0;
113 int parallel_ref_levels = 0;
114 bool visualization = true;
115 bool visit = true;
116
117 Array<int> dbcs;
118 Array<int> nbcs;
119
120 Vector dbcv;
121 Vector nbcv;
122
123 bool dbcg = false;
124
125 OptionsParser args(argc, argv);
126 args.AddOption(&mesh_file, "-m", "--mesh",
127 "Mesh file to use.");
128 args.AddOption(&order, "-o", "--order",
129 "Finite element order (polynomial degree).");
130 args.AddOption(&serial_ref_levels, "-rs", "--serial-ref-levels",
131 "Number of serial refinement levels.");
132 args.AddOption(&parallel_ref_levels, "-rp", "--parallel-ref-levels",
133 "Number of parallel refinement levels.");
134 args.AddOption(&e_uniform_, "-uebc", "--uniform-e-bc",
135 "Specify if the three components of the constant "
136 "electric field");
137 args.AddOption(&pw_eps_, "-pwe", "--piecewise-eps",
138 "Piecewise values of Permittivity");
139 args.AddOption(&ds_params_, "-ds", "--dielectric-sphere-params",
140 "Center, Radius, and Permittivity of Dielectric Sphere");
141 args.AddOption(&cs_params_, "-cs", "--charged-sphere-params",
142 "Center, Radius, and Total Charge of Charged Sphere");
143 args.AddOption(&pc_params_, "-pc", "--point-charge-params",
144 "Charges and locations of Point Charges");
145 args.AddOption(&vp_params_, "-vp", "--voltaic-pile-params",
146 "Axis End Points, Radius, and "
147 "Polarization of Cylindrical Voltaic Pile");
148 args.AddOption(&dbcs, "-dbcs", "--dirichlet-bc-surf",
149 "Dirichlet Boundary Condition Surfaces");
150 args.AddOption(&dbcv, "-dbcv", "--dirichlet-bc-vals",
151 "Dirichlet Boundary Condition Values");
154 "Dirichlet Boundary Condition Gradient (phi = -z)");
155 args.AddOption(&nbcs, "-nbcs", "--neumann-bc-surf",
156 "Neumann Boundary Condition Surfaces");
157 args.AddOption(&nbcv, "-nbcv", "--neumann-bc-vals",
158 "Neumann Boundary Condition Values");
159 args.AddOption(&maxit, "-maxit", "--max-amr-iterations",
160 "Max number of iterations in the main AMR loop.");
161 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
162 "--no-visualization",
163 "Enable or disable GLVis visualization.");
164 args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit",
165 "Enable or disable VisIt visualization.");
166 args.Parse();
167 if (!args.Good())
168 {
169 if (Mpi::Root())
170 {
171 args.PrintUsage(cout);
172 }
173 return 1;
174 }
175 if (Mpi::Root())
176 {
177 args.PrintOptions(cout);
178 }
179
180 // Read the (serial) mesh from the given mesh file on all processors. We
181 // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
182 // and volume meshes with the same code.
183 Mesh *mesh = new Mesh(mesh_file, 1, 1);
184 int sdim = mesh->SpaceDimension();
185
186 if (Mpi::Root())
187 {
188 cout << "Starting initialization." << endl;
189 }
190
191 // Project a NURBS mesh to a piecewise-quadratic curved mesh
192 if (mesh->NURBSext)
193 {
194 mesh->UniformRefinement();
195 if (serial_ref_levels > 0) { serial_ref_levels--; }
196
197 mesh->SetCurvature(2);
198 }
199
200 // Ensure that quad and hex meshes are treated as non-conforming.
201 mesh->EnsureNCMesh();
202
203 // Refine the serial mesh on all processors to increase the resolution. In
204 // this example we do 'ref_levels' of uniform refinement. NURBS meshes are
205 // refined at least twice, as they are typically coarse.
206 for (int l = 0; l < serial_ref_levels; l++)
207 {
208 mesh->UniformRefinement();
209 }
210
211 // Define a parallel mesh by a partitioning of the serial mesh. Refine
212 // this mesh further in parallel to increase the resolution. Once the
213 // parallel mesh is defined, the serial mesh can be deleted.
214 ParMesh pmesh(MPI_COMM_WORLD, *mesh);
215 delete mesh;
216
217 // Refine this mesh in parallel to increase the resolution.
218 int par_ref_levels = parallel_ref_levels;
219 for (int l = 0; l < par_ref_levels; l++)
220 {
221 pmesh.UniformRefinement();
222 }
223 // Make sure tet-only meshes are marked for local refinement.
224 pmesh.Finalize(true);
225
226 // If the gradient bc was selected but the E field was not specified
227 // set a default vector value.
228 if ( dbcg && e_uniform_.Size() != sdim )
229 {
230 e_uniform_.SetSize(sdim);
231 e_uniform_ = 0.0;
232 e_uniform_(sdim-1) = 1.0;
233 }
234
235 // If values for Dirichlet BCs were not set assume they are zero
236 if ( dbcv.Size() < dbcs.Size() && !dbcg )
237 {
238 dbcv.SetSize(dbcs.Size());
239 dbcv = 0.0;
240 }
241
242 // If values for Neumann BCs were not set assume they are zero
243 if ( nbcv.Size() < nbcs.Size() )
244 {
245 nbcv.SetSize(nbcs.Size());
246 nbcv = 0.0;
247 }
248
249 // Create a coefficient describing the dielectric permittivity
251
252 // Create the Electrostatic solver
253 VoltaSolver Volta(pmesh, order, dbcs, dbcv, nbcs, nbcv, *epsCoef,
254 ( e_uniform_.Size() > 0 ) ? phi_bc_uniform : NULL,
255 ( cs_params_.Size() > 0 ) ? charged_sphere : NULL,
256 ( vp_params_.Size() > 0 ) ? voltaic_pile : NULL,
257 pc_params_);
258
259 // Initialize GLVis visualization
260 if (visualization)
261 {
262 Volta.InitializeGLVis();
263 }
264
265 // Initialize VisIt visualization
266 VisItDataCollection visit_dc("Volta-AMR-Parallel", &pmesh);
267
268 if ( visit )
269 {
270 Volta.RegisterVisItFields(visit_dc);
271 }
272 if (Mpi::Root()) { cout << "Initialization done." << endl; }
273
274 // The main AMR loop. In each iteration we solve the problem on the current
275 // mesh, visualize the solution, estimate the error on all elements, refine
276 // the worst elements and update all objects to work with the new mesh. We
277 // refine until the maximum number of dofs in the nodal finite element space
278 // reaches 10 million.
279 const int max_dofs = 10000000;
280 for (int it = 1; it <= maxit; it++)
281 {
282 if (Mpi::Root())
283 {
284 cout << "\nAMR Iteration " << it << endl;
285 }
286
287 // Display the current number of DoFs in each finite element space
288 Volta.PrintSizes();
289
290 // Assemble all forms
291 Volta.Assemble();
292
293 // Solve the system and compute any auxiliary fields
294 Volta.Solve();
295
296 // Determine the current size of the linear system
297 int prob_size = Volta.GetProblemSize();
298
299 // Write fields to disk for VisIt
300 if ( visit )
301 {
302 Volta.WriteVisItFields(it);
303 }
304
305 // Send the solution by socket to a GLVis server.
306 if (visualization)
307 {
308 Volta.DisplayToGLVis();
309 }
310
311 if (Mpi::Root())
312 {
313 cout << "AMR iteration " << it << " complete." << endl;
314 }
315
316 // Check stopping criteria
317 if (prob_size > max_dofs)
318 {
319 if (Mpi::Root())
320 {
321 cout << "Reached maximum number of dofs, exiting..." << endl;
322 }
323 break;
324 }
325 if (it == maxit)
326 {
327 break;
328 }
329
330 // Wait for user input. Ask every 10th iteration.
331 char c = 'c';
332 if (Mpi::Root() && (it % 10 == 0))
333 {
334 cout << "press (q)uit or (c)ontinue --> " << flush;
335 cin >> c;
336 }
337 MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
338
339 if (c != 'c')
340 {
341 break;
342 }
343
344 // Estimate element errors using the Zienkiewicz-Zhu error estimator.
345 Vector errors(pmesh.GetNE());
346 Volta.GetErrorEstimates(errors);
347
348 real_t local_max_err = errors.Max();
349 real_t global_max_err;
350 MPI_Allreduce(&local_max_err, &global_max_err, 1,
351 MPITypeMap<real_t>::mpi_type, MPI_MAX, pmesh.GetComm());
352
353 // Refine the elements whose error is larger than a fraction of the
354 // maximum element error.
355 const real_t frac = 0.7;
356 real_t threshold = frac * global_max_err;
357 if (Mpi::Root()) { cout << "Refining ..." << endl; }
358 pmesh.RefineByError(errors, threshold);
359
360 // Update the electrostatic solver to reflect the new state of the mesh.
361 Volta.Update();
362
363 if (pmesh.Nonconforming() && Mpi::WorldSize() > 1)
364 {
365 if (Mpi::Root()) { cout << "Rebalancing ..." << endl; }
366 pmesh.Rebalance();
367
368 // Update again after rebalancing
369 Volta.Update();
370 }
371 }
372
373 delete epsCoef;
374
375 return 0;
376}
377
378// Print the Volta ascii logo to the given ostream
379void display_banner(ostream & os)
380{
381 os << " ____ ____ __ __ " << endl
382 << " \\ \\ / /___ | |_/ |______ " << endl
383 << " \\ Y / _ \\| |\\ __\\__ \\ " << endl
384 << " \\ ( <_> ) |_| | / __ \\_ " << endl
385 << " \\___/ \\____/|____/__| (____ / " << endl
386 << " \\/ " << endl << flush;
387}
388
389// The Permittivity is a required coefficient which may be defined in
390// various ways so we'll determine the appropriate coefficient type here.
393{
394 Coefficient * coef = NULL;
395
396 if ( ds_params_.Size() > 0 )
397 {
399 }
400 else if ( pw_eps_.Size() > 0 )
401 {
402 coef = new PWConstCoefficient(pw_eps_);
403 }
404 else
405 {
406 coef = new ConstantCoefficient(epsilon0_);
407 }
408
409 return coef;
410}
411
412// A sphere with constant permittivity. The sphere has a radius,
413// center, and permittivity specified on the command line and stored
414// in ds_params_.
416{
417 real_t r2 = 0.0;
418
419 for (int i=0; i<x.Size(); i++)
420 {
421 r2 += (x(i)-ds_params_(i))*(x(i)-ds_params_(i));
422 }
423
424 if ( sqrt(r2) <= ds_params_(x.Size()) )
425 {
426 return ds_params_(x.Size()+1) * epsilon0_;
427 }
428 return epsilon0_;
429}
430
431// A sphere with constant charge density. The sphere has a radius,
432// center, and total charge specified on the command line and stored
433// in cs_params_.
435{
436 real_t r2 = 0.0;
437 real_t rho = 0.0;
438
439 if ( cs_params_(x.Size()) > 0.0 )
440 {
441 switch ( x.Size() )
442 {
443 case 2:
444 rho = cs_params_(x.Size()+1) /
445 (M_PI * pow(cs_params_(x.Size()), 2));
446 break;
447 case 3:
448 rho = 0.75 * cs_params_(x.Size()+1) /
449 (M_PI * pow(cs_params_(x.Size()), 3));
450 break;
451 default:
452 rho = 0.0;
453 }
454 }
455
456 for (int i=0; i<x.Size(); i++)
457 {
458 r2 += (x(i) - cs_params_(i)) * (x(i) - cs_params_(i));
459 }
460
461 if ( sqrt(r2) <= cs_params_(x.Size()) )
462 {
463 return rho;
464 }
465 return 0.0;
466}
467
468// A Cylindrical Rod of constant polarization. The cylinder has two
469// axis end points, a radius, and a constant electric polarization oriented
470// along the axis.
471void voltaic_pile(const Vector &x, Vector &p)
472{
473 p.SetSize(x.Size());
474 p = 0.0;
475
476 Vector a(x.Size()); // Normalized Axis vector
477 Vector xu(x.Size()); // x vector relative to the axis end-point
478
479 xu = x;
480
481 for (int i=0; i<x.Size(); i++)
482 {
483 xu[i] -= vp_params_[i];
484 a[i] = vp_params_[x.Size()+i] - vp_params_[i];
485 }
486
487 real_t h = a.Norml2();
488
489 if ( h == 0.0 )
490 {
491 return;
492 }
493
494 real_t r = vp_params_[2 * x.Size()];
495 real_t xa = xu * a;
496
497 if ( h > 0.0 )
498 {
499 xu.Add(-xa / (h * h), a);
500 }
501
502 real_t xp = xu.Norml2();
503
504 if ( xa >= 0.0 && xa <= h*h && xp <= r )
505 {
506 p.Add(vp_params_[2 * x.Size() + 1] / h, a);
507 }
508}
509
510// To produce a uniform electric field the potential can be set
511// to (- Ex x - Ey y - Ez z).
513{
514 real_t phi = 0.0;
515
516 for (int i=0; i<x.Size(); i++)
517 {
518 phi -= x(i) * e_uniform_(i);
519 }
520
521 return phi;
522}
int Size() const
Return the logical size of the array.
Definition array.hpp:144
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
A coefficient that is constant across space and time.
A general function coefficient.
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
bool Nonconforming() const
Definition mesh.hpp:2229
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1163
bool RefineByError(const Array< real_t > &elem_error, real_t threshold, int nonconforming=-1, int nc_limit=0)
Definition mesh.cpp:10695
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition mesh.cpp:10626
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 int WorldSize()
Return the size of MPI_COMM_WORLD.
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.
A piecewise constant coefficient with the constants keyed off the element attribute numbers.
Class for parallel meshes.
Definition pmesh.hpp:34
MPI_Comm GetComm() const
Definition pmesh.hpp:402
void Rebalance()
Definition pmesh.cpp:4009
void Finalize(bool refine=false, bool fix_orientation=false) override
Finalize the construction of a general Mesh.
Definition pmesh.cpp:1524
Vector data type.
Definition vector.hpp:80
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:832
real_t Max() const
Returns the maximal element of the vector.
Definition vector.cpp:922
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 RegisterVisItFields(VisItDataCollection &visit_dc)
int main()
real_t a
Definition lissajous.cpp:41
float real_t
Definition config.hpp:43
real_t p(const Vector &x, real_t t)
Helper struct to convert a C++ type to an MPI type.
Coefficient * SetupPermittivityCoefficient()
Definition volta.cpp:392
real_t phi_bc_uniform(const Vector &)
Definition volta.cpp:512
real_t dielectric_sphere(const Vector &)
Definition volta.cpp:415
real_t charged_sphere(const Vector &)
Definition volta.cpp:434
void display_banner(ostream &os)
Definition volta.cpp:379
void voltaic_pile(const Vector &, Vector &)
Definition volta.cpp:471