MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
volta.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 
69 using namespace std;
70 using namespace mfem;
71 using namespace mfem::electromagnetics;
72 
73 // Permittivity Functions
75 
76 static Vector pw_eps_(0); // Piecewise permittivity values
77 static Vector ds_params_(0); // Center, Radius, and Permittivity
78 // of dielectric sphere
79 double dielectric_sphere(const Vector &);
80 
81 // Charge Density Function
82 static Vector cs_params_(0); // Center, Radius, and Total Charge
83 // of charged sphere
84 double charged_sphere(const Vector &);
85 
86 // Point Charges
87 static Vector pc_params_(0); // Point charge locations and charges
88 
89 // Polarization
90 static Vector vp_params_(0); // Axis Start, Axis End, Cylinder Radius,
91 // and Polarization Magnitude
92 void voltaic_pile(const Vector &, Vector &);
93 
94 // Phi Boundary Condition
95 static Vector e_uniform_(0);
96 double phi_bc_uniform(const Vector &);
97 
98 // Prints the program's logo to the given output stream
99 void display_banner(ostream & os);
100 
101 int 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");
152  args.AddOption(&dbcg, "-dbcg", "--dirichlet-bc-gradient",
153  "-no-dbcg", "--no-dirichlet-bc-gradient",
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  double local_max_err = errors.Max();
349  double global_max_err;
350  MPI_Allreduce(&local_max_err, &global_max_err, 1,
351  MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
352 
353  // Refine the elements whose error is larger than a fraction of the
354  // maximum element error.
355  const double frac = 0.7;
356  double 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
379 void 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.
391 Coefficient *
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_.
415 double dielectric_sphere(const Vector &x)
416 {
417  double 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_.
434 double charged_sphere(const Vector &x)
435 {
436  double r2 = 0.0;
437  double 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.
471 void 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  double h = a.Norml2();
488 
489  if ( h == 0.0 )
490  {
491  return;
492  }
493 
494  double r = vp_params_[2 * x.Size()];
495  double xa = xu * a;
496 
497  if ( h > 0.0 )
498  {
499  xu.Add(-xa / (h * h), a);
500  }
501 
502  double 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).
512 double phi_bc_uniform(const Vector &x)
513 {
514  double 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 }
double phi_bc_uniform(const Vector &)
Definition: volta.cpp:512
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
double dielectric_sphere(const Vector &)
Definition: maxwell.cpp:382
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:923
void RegisterVisItFields(VisItDataCollection &visit_dc)
void voltaic_pile(const Vector &, Vector &)
Definition: volta.cpp:471
void GetErrorEstimates(Vector &errors)
void Rebalance()
Definition: pmesh.cpp:3989
bool Nonconforming() const
Definition: mesh.hpp:1651
void Finalize(bool refine=false, bool fix_orientation=false) override
Finalize the construction of a general Mesh.
Definition: pmesh.cpp:1516
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition: mesh.cpp:9473
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9816
Data collection with VisIt I/O routines.
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:5343
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
double charged_sphere(const Vector &)
Definition: volta.cpp:434
int SpaceDimension() const
Definition: mesh.hpp:1007
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
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
bool RefineByError(const Array< double > &elem_error, double threshold, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:9542
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:272
A piecewise constant coefficient with the constants keyed off the element attribute numbers...
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:252
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
Coefficient * SetupPermittivityCoefficient()
Definition: volta.cpp:392
void display_banner(ostream &os)
Definition: joule.cpp:778
Class for parallel meshes.
Definition: pmesh.hpp:32
int main()
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150