MFEM  v4.3.0 Finite element discretization library
volta.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2021, 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 // 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_Session mpi(argc, argv);
104
105  if ( mpi.Root() ) { display_banner(cout); }
106
107  // Parse command-line options.
108  const char *mesh_file = "../../data/ball-nurbs.mesh";
109  int order = 1;
110  int maxit = 100;
111  int serial_ref_levels = 0;
112  int parallel_ref_levels = 0;
113  bool visualization = true;
114  bool visit = true;
115
116  Array<int> dbcs;
117  Array<int> nbcs;
118
119  Vector dbcv;
120  Vector nbcv;
121
122  bool dbcg = false;
123
124  OptionsParser args(argc, argv);
126  "Mesh file to use.");
128  "Finite element order (polynomial degree).");
130  "Number of serial refinement levels.");
132  "Number of parallel refinement levels.");
134  "Specify if the three components of the constant "
135  "electric field");
137  "Piecewise values of Permittivity");
139  "Center, Radius, and Permittivity of Dielectric Sphere");
141  "Center, Radius, and Total Charge of Charged Sphere");
143  "Charges and locations of Point Charges");
145  "Axis End Points, Radius, and "
146  "Polarization of Cylindrical Voltaic Pile");
148  "Dirichlet Boundary Condition Surfaces");
150  "Dirichlet Boundary Condition Values");
153  "Dirichlet Boundary Condition Gradient (phi = -z)");
155  "Neumann Boundary Condition Surfaces");
157  "Neumann Boundary Condition Values");
159  "Max number of iterations in the main AMR loop.");
161  "--no-visualization",
162  "Enable or disable GLVis visualization.");
163  args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit",
164  "Enable or disable VisIt visualization.");
165  args.Parse();
166  if (!args.Good())
167  {
168  if (mpi.Root())
169  {
170  args.PrintUsage(cout);
171  }
172  return 1;
173  }
174  if (mpi.Root())
175  {
176  args.PrintOptions(cout);
177  }
178
179  // Read the (serial) mesh from the given mesh file on all processors. We
180  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
181  // and volume meshes with the same code.
182  Mesh *mesh = new Mesh(mesh_file, 1, 1);
183  int sdim = mesh->SpaceDimension();
184
185  if (mpi.Root())
186  {
187  cout << "Starting initialization." << endl;
188  }
189
190  // Project a NURBS mesh to a piecewise-quadratic curved mesh
191  if (mesh->NURBSext)
192  {
193  mesh->UniformRefinement();
194  if (serial_ref_levels > 0) { serial_ref_levels--; }
195
196  mesh->SetCurvature(2);
197  }
198
199  // Ensure that quad and hex meshes are treated as non-conforming.
200  mesh->EnsureNCMesh();
201
202  // Refine the serial mesh on all processors to increase the resolution. In
203  // this example we do 'ref_levels' of uniform refinement. NURBS meshes are
204  // refined at least twice, as they are typically coarse.
205  for (int l = 0; l < serial_ref_levels; l++)
206  {
207  mesh->UniformRefinement();
208  }
209
210  // Define a parallel mesh by a partitioning of the serial mesh. Refine
211  // this mesh further in parallel to increase the resolution. Once the
212  // parallel mesh is defined, the serial mesh can be deleted.
213  ParMesh pmesh(MPI_COMM_WORLD, *mesh);
214  delete mesh;
215
216  // Refine this mesh in parallel to increase the resolution.
217  int par_ref_levels = parallel_ref_levels;
218  for (int l = 0; l < par_ref_levels; l++)
219  {
220  pmesh.UniformRefinement();
221  }
222  // Make sure tet-only meshes are marked for local refinement.
223  pmesh.Finalize(true);
224
225  // If the gradient bc was selected but the E field was not specified
226  // set a default vector value.
227  if ( dbcg && e_uniform_.Size() != sdim )
228  {
229  e_uniform_.SetSize(sdim);
230  e_uniform_ = 0.0;
231  e_uniform_(sdim-1) = 1.0;
232  }
233
234  // If values for Dirichlet BCs were not set assume they are zero
235  if ( dbcv.Size() < dbcs.Size() && !dbcg )
236  {
237  dbcv.SetSize(dbcs.Size());
238  dbcv = 0.0;
239  }
240
241  // If values for Neumann BCs were not set assume they are zero
242  if ( nbcv.Size() < nbcs.Size() )
243  {
244  nbcv.SetSize(nbcs.Size());
245  nbcv = 0.0;
246  }
247
248  // Create a coefficient describing the dielectric permittivity
250
251  // Create the Electrostatic solver
252  VoltaSolver Volta(pmesh, order, dbcs, dbcv, nbcs, nbcv, *epsCoef,
253  ( e_uniform_.Size() > 0 ) ? phi_bc_uniform : NULL,
254  ( cs_params_.Size() > 0 ) ? charged_sphere : NULL,
255  ( vp_params_.Size() > 0 ) ? voltaic_pile : NULL,
256  pc_params_);
257
258  // Initialize GLVis visualization
259  if (visualization)
260  {
261  Volta.InitializeGLVis();
262  }
263
264  // Initialize VisIt visualization
265  VisItDataCollection visit_dc("Volta-AMR-Parallel", &pmesh);
266
267  if ( visit )
268  {
269  Volta.RegisterVisItFields(visit_dc);
270  }
271  if (mpi.Root()) { cout << "Initialization done." << endl; }
272
273  // The main AMR loop. In each iteration we solve the problem on the current
274  // mesh, visualize the solution, estimate the error on all elements, refine
275  // the worst elements and update all objects to work with the new mesh. We
276  // refine until the maximum number of dofs in the nodal finite element space
277  // reaches 10 million.
278  const int max_dofs = 10000000;
279  for (int it = 1; it <= maxit; it++)
280  {
281  if (mpi.Root())
282  {
283  cout << "\nAMR Iteration " << it << endl;
284  }
285
286  // Display the current number of DoFs in each finite element space
287  Volta.PrintSizes();
288
289  // Assemble all forms
290  Volta.Assemble();
291
292  // Solve the system and compute any auxiliary fields
293  Volta.Solve();
294
295  // Determine the current size of the linear system
296  int prob_size = Volta.GetProblemSize();
297
298  // Write fields to disk for VisIt
299  if ( visit )
300  {
301  Volta.WriteVisItFields(it);
302  }
303
304  // Send the solution by socket to a GLVis server.
305  if (visualization)
306  {
307  Volta.DisplayToGLVis();
308  }
309
310  if (mpi.Root())
311  {
312  cout << "AMR iteration " << it << " complete." << endl;
313  }
314
315  // Check stopping criteria
316  if (prob_size > max_dofs)
317  {
318  if (mpi.Root())
319  {
320  cout << "Reached maximum number of dofs, exiting..." << endl;
321  }
322  break;
323  }
324  if (it == maxit)
325  {
326  break;
327  }
328
329  // Wait for user input. Ask every 10th iteration.
330  char c = 'c';
331  if (mpi.Root() && (it % 10 == 0))
332  {
333  cout << "press (q)uit or (c)ontinue --> " << flush;
334  cin >> c;
335  }
336  MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
337
338  if (c != 'c')
339  {
340  break;
341  }
342
343  // Estimate element errors using the Zienkiewicz-Zhu error estimator.
344  Vector errors(pmesh.GetNE());
345  Volta.GetErrorEstimates(errors);
346
347  double local_max_err = errors.Max();
348  double global_max_err;
349  MPI_Allreduce(&local_max_err, &global_max_err, 1,
350  MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
351
352  // Refine the elements whose error is larger than a fraction of the
353  // maximum element error.
354  const double frac = 0.7;
355  double threshold = frac * global_max_err;
356  if (mpi.Root()) { cout << "Refining ..." << endl; }
357  pmesh.RefineByError(errors, threshold);
358
359  // Update the electrostatic solver to reflect the new state of the mesh.
360  Volta.Update();
361
362  if (pmesh.Nonconforming() && mpi.WorldSize() > 1)
363  {
364  if (mpi.Root()) { cout << "Rebalancing ..." << endl; }
365  pmesh.Rebalance();
366
367  // Update again after rebalancing
368  Volta.Update();
369  }
370  }
371
372  delete epsCoef;
373
374  return 0;
375 }
376
377 // Print the Volta ascii logo to the given ostream
378 void display_banner(ostream & os)
379 {
380  os << " ____ ____ __ __ " << endl
381  << " \\ \\ / /___ | |_/ |______ " << endl
382  << " \\ Y / _ \\| |\\ __\\__ \\ " << endl
383  << " \\ ( <_> ) |_| | / __ \\_ " << endl
384  << " \\___/ \\____/|____/__| (____ / " << endl
385  << " \\/ " << endl << flush;
386 }
387
388 // The Permittivity is a required coefficient which may be defined in
389 // various ways so we'll determine the appropriate coefficient type here.
390 Coefficient *
392 {
393  Coefficient * coef = NULL;
394
395  if ( ds_params_.Size() > 0 )
396  {
398  }
399  else if ( pw_eps_.Size() > 0 )
400  {
401  coef = new PWConstCoefficient(pw_eps_);
402  }
403  else
404  {
405  coef = new ConstantCoefficient(epsilon0_);
406  }
407
408  return coef;
409 }
410
411 // A sphere with constant permittivity. The sphere has a radius,
412 // center, and permittivity specified on the command line and stored
413 // in ds_params_.
414 double dielectric_sphere(const Vector &x)
415 {
416  double r2 = 0.0;
417
418  for (int i=0; i<x.Size(); i++)
419  {
420  r2 += (x(i)-ds_params_(i))*(x(i)-ds_params_(i));
421  }
422
423  if ( sqrt(r2) <= ds_params_(x.Size()) )
424  {
425  return ds_params_(x.Size()+1) * epsilon0_;
426  }
427  return epsilon0_;
428 }
429
430 // A sphere with constant charge density. The sphere has a radius,
431 // center, and total charge specified on the command line and stored
432 // in cs_params_.
433 double charged_sphere(const Vector &x)
434 {
435  double r2 = 0.0;
436  double rho = 0.0;
437
438  if ( cs_params_(x.Size()) > 0.0 )
439  {
440  switch ( x.Size() )
441  {
442  case 2:
443  rho = cs_params_(x.Size()+1) /
444  (M_PI * pow(cs_params_(x.Size()), 2));
445  break;
446  case 3:
447  rho = 0.75 * cs_params_(x.Size()+1) /
448  (M_PI * pow(cs_params_(x.Size()), 3));
449  break;
450  default:
451  rho = 0.0;
452  }
453  }
454
455  for (int i=0; i<x.Size(); i++)
456  {
457  r2 += (x(i) - cs_params_(i)) * (x(i) - cs_params_(i));
458  }
459
460  if ( sqrt(r2) <= cs_params_(x.Size()) )
461  {
462  return rho;
463  }
464  return 0.0;
465 }
466
467 // A Cylindrical Rod of constant polarization. The cylinder has two
468 // axis end points, a radius, and a constant electric polarization oriented
469 // along the axis.
470 void voltaic_pile(const Vector &x, Vector &p)
471 {
472  p.SetSize(x.Size());
473  p = 0.0;
474
475  Vector a(x.Size()); // Normalized Axis vector
476  Vector xu(x.Size()); // x vector relative to the axis end-point
477
478  xu = x;
479
480  for (int i=0; i<x.Size(); i++)
481  {
482  xu[i] -= vp_params_[i];
483  a[i] = vp_params_[x.Size()+i] - vp_params_[i];
484  }
485
486  double h = a.Norml2();
487
488  if ( h == 0.0 )
489  {
490  return;
491  }
492
493  double r = vp_params_[2 * x.Size()];
494  double xa = xu * a;
495
496  if ( h > 0.0 )
497  {
498  xu.Add(-xa / (h * h), a);
499  }
500
501  double xp = xu.Norml2();
502
503  if ( xa >= 0.0 && xa <= h*h && xp <= r )
504  {
505  p.Add(vp_params_[2 * x.Size() + 1] / h, a);
506  }
507 }
508
509 // To produce a uniform electric field the potential can be set
510 // to (- Ex x - Ey y - Ez z).
511 double phi_bc_uniform(const Vector &x)
512 {
513  double phi = 0.0;
514
515  for (int i=0; i<x.Size(); i++)
516  {
517  phi -= x(i) * e_uniform_(i);
518  }
519
520  return phi;
521 }
double phi_bc_uniform(const Vector &)
Definition: volta.cpp:511
int WorldSize() const
Return MPI_COMM_WORLD&#39;s size.
int Size() const
Return the logical size of the array.
Definition: array.hpp:134
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
double dielectric_sphere(const Vector &)
Definition: maxwell.cpp:381
int Size() const
Returns the size of the vector.
Definition: vector.hpp:190
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:846
void RegisterVisItFields(VisItDataCollection &visit_dc)
void voltaic_pile(const Vector &, Vector &)
Definition: volta.cpp:470
void GetErrorEstimates(Vector &errors)
void Rebalance()
Definition: pmesh.cpp:3602
bool Nonconforming() const
Definition: mesh.hpp:1367
A simple convenience class that calls MPI_Init() at construction and MPI_Finalize() at destruction...
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:150
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition: mesh.cpp:8800
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
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:4882
bool Root() const
Return true if WorldRank() == 0.
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:457
double charged_sphere(const Vector &)
Definition: volta.cpp:433
int SpaceDimension() const
Definition: mesh.hpp:912
MPI_Comm GetComm() const
Definition: pmesh.hpp:276
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:39
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:8869
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:206
A piecewise constant coefficient with the constants keyed off the element attribute numbers...
Definition: coefficient.hpp:94
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:243
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:327
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
Definition: pmesh.cpp:1510
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
Coefficient * SetupPermittivityCoefficient()
Definition: volta.cpp:391
void display_banner(ostream &os)
Definition: joule.cpp:780
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