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