MFEM  v3.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, 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 // By default the sources and fields are all zero:
54 // mpirun -np 4 volta
55 
56 #include "volta_solver.hpp"
57 #include <fstream>
58 #include <iostream>
59 
60 using namespace std;
61 using namespace mfem;
62 using namespace mfem::electromagnetics;
63 
64 // Permittivity Function
65 static Vector ds_params_(0); // Center, Radius, and Permittivity
66 // of dielectric sphere
67 double dielectric_sphere(const Vector &);
68 
69 // Charge Density Function
70 static Vector cs_params_(0); // Center, Radius, and Total Charge
71 // of charged sphere
72 double charged_sphere(const Vector &);
73 
74 // Polarization
75 static Vector vp_params_(0); // Axis Start, Axis End, Cylinder Radius,
76 // and Polarization Magnitude
77 void voltaic_pile(const Vector &, Vector &);
78 
79 // Phi Boundary Condition
80 static Vector e_uniform_(0);
81 double phi_bc_uniform(const Vector &);
82 
83 // Prints the program's logo to the given output stream
84 void display_banner(ostream & os);
85 
86 int main(int argc, char *argv[])
87 {
88  // Initialize MPI.
89  int num_procs, myid;
90  MPI_Init(&argc, &argv);
91  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
92  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
93 
94  if ( myid == 0 ) { display_banner(cout); }
95 
96  // Parse command-line options.
97  const char *mesh_file = "butterfly_3d.mesh";
98  int order = 1;
99  int sr = 0, pr = 0;
100  bool visualization = true;
101  bool visit = true;
102 
103  Array<int> dbcs;
104  Array<int> nbcs;
105 
106  Vector dbcv;
107  Vector nbcv;
108 
109  bool dbcg = false;
110 
111  OptionsParser args(argc, argv);
112  args.AddOption(&mesh_file, "-m", "--mesh",
113  "Mesh file to use.");
114  args.AddOption(&order, "-o", "--order",
115  "Finite element order (polynomial degree).");
116  args.AddOption(&sr, "-rs", "--serial-ref-levels",
117  "Number of serial refinement levels.");
118  args.AddOption(&pr, "-rp", "--parallel-ref-levels",
119  "Number of parallel refinement levels.");
120  args.AddOption(&e_uniform_, "-uebc", "--uniform-e-bc",
121  "Specify if the three components of the constant electric field");
122  args.AddOption(&ds_params_, "-ds", "--dielectric-sphere-params",
123  "Center, Radius, and Permittivity of Dielectric Sphere");
124  args.AddOption(&cs_params_, "-cs", "--charged-sphere-params",
125  "Center, Radius, and Total Charge of Charged Sphere");
126  args.AddOption(&vp_params_, "-vp", "--voltaic-pile-params",
127  "Axis End Points, Radius, and Polarization of Cylindrical Voltaic Pile");
128  args.AddOption(&dbcs, "-dbcs", "--dirichlet-bc-surf",
129  "Dirichlet Boundary Condition Surfaces");
130  args.AddOption(&dbcv, "-dbcv", "--dirichlet-bc-vals",
131  "Dirichlet Boundary Condition Values");
132  args.AddOption(&dbcg, "-dbcg", "--dirichlet-bc-gradient",
133  "-no-dbcg", "--no-dirichlet-bc-gradient",
134  "Dirichlet Boundary Condition Gradient (phi = -z)");
135  args.AddOption(&nbcs, "-nbcs", "--neumann-bc-surf",
136  "Neumann Boundary Condition Surfaces");
137  args.AddOption(&nbcv, "-nbcv", "--neumann-bc-vals",
138  "Neumann Boundary Condition Values");
139  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
140  "--no-visualization",
141  "Enable or disable GLVis visualization.");
142  args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit",
143  "Enable or disable VisIt visualization.");
144  args.Parse();
145  if (!args.Good())
146  {
147  if (myid == 0)
148  {
149  args.PrintUsage(cout);
150  }
151  MPI_Finalize();
152  return 1;
153  }
154  if (myid == 0)
155  {
156  args.PrintOptions(cout);
157  }
158 
159  // Read the (serial) mesh from the given mesh file on all processors. We
160  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
161  // and volume meshes with the same code.
162  Mesh *mesh;
163  ifstream imesh(mesh_file);
164  if (!imesh)
165  {
166  if (myid == 0)
167  {
168  cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
169  }
170  MPI_Finalize();
171  return 2;
172  }
173  mesh = new Mesh(imesh, 1, 1);
174  imesh.close();
175 
176  int sdim = mesh->SpaceDimension();
177 
178  // Refine the serial mesh on all processors to increase the resolution. In
179  // this example we do 'ref_levels' of uniform refinement. NURBS meshes are
180  // refined at least twice, as they are typically coarse.
181  if (myid == 0) { cout << "Starting initialization." << endl; }
182  {
183  int ref_levels = sr;
184  if (mesh->NURBSext && ref_levels < 2)
185  {
186  ref_levels = 2;
187  }
188  for (int l = 0; l < ref_levels; l++)
189  {
190  mesh->UniformRefinement();
191  }
192  }
193 
194  // Project a NURBS mesh to a piecewise-quadratic curved mesh. Make sure that
195  // the mesh is non-conforming.
196  if (mesh->NURBSext)
197  {
198  mesh->SetCurvature(2);
199  }
200  mesh->EnsureNCMesh();
201 
202  // Define a parallel mesh by a partitioning of the serial mesh. Refine
203  // this mesh further in parallel to increase the resolution. Once the
204  // parallel mesh is defined, the serial mesh can be deleted.
205  ParMesh pmesh(MPI_COMM_WORLD, *mesh);
206  delete mesh;
207 
208  // Refine this mesh in parallel to increase the resolution.
209  int par_ref_levels = pr;
210  for (int l = 0; l < par_ref_levels; l++)
211  {
212  pmesh.UniformRefinement();
213  }
214 
215  // If the gradient bc was selected but the E field was not specified
216  // set a default vector value.
217  if ( dbcg && e_uniform_.Size() != sdim )
218  {
219  e_uniform_.SetSize(sdim);
220  e_uniform_ = 0.0;
221  e_uniform_(sdim-1) = 1.0;
222  }
223 
224  // If values for Dirichlet BCs were not set assume they are zero
225  if ( dbcv.Size() < dbcs.Size() && !dbcg )
226  {
227  dbcv.SetSize(dbcs.Size());
228  dbcv = 0.0;
229  }
230 
231  // If values for Neumann BCs were not set assume they are zero
232  if ( nbcv.Size() < nbcs.Size() )
233  {
234  nbcv.SetSize(nbcs.Size());
235  nbcv = 0.0;
236  }
237 
238  // Create the Electrostatic solver
239  VoltaSolver Volta(pmesh, order, dbcs, dbcv, nbcs, nbcv,
240  ( ds_params_.Size() > 0 ) ? dielectric_sphere : NULL,
241  ( e_uniform_.Size() > 0 ) ? phi_bc_uniform : NULL,
242  ( cs_params_.Size() > 0 ) ? charged_sphere : NULL,
243  ( vp_params_.Size() > 0 ) ? voltaic_pile : NULL);
244 
245  // Initialize GLVis visualization
246  if (visualization)
247  {
248  Volta.InitializeGLVis();
249  }
250 
251  // Initialize VisIt visualization
252  VisItDataCollection visit_dc("Volta-AMR-Parallel", &pmesh);
253 
254  if ( visit )
255  {
256  Volta.RegisterVisItFields(visit_dc);
257  }
258  if (myid == 0) { cout << "Initialization done." << endl; }
259 
260  // The main AMR loop. In each iteration we solve the problem on the current
261  // mesh, visualize the solution, estimate the error on all elements, refine
262  // the worst elements and update all objects to work with the new mesh. We
263  // refine until the maximum number of dofs in the nodal finite element space
264  // reaches 10 million.
265  const int max_dofs = 10000000;
266  for (int it = 1; it <= 100; it++)
267  {
268  if (myid == 0)
269  {
270  cout << "\nAMR Iteration " << it << endl;
271  }
272 
273  // Display the current number of DoFs in each finite element space
274  Volta.PrintSizes();
275 
276  // Solve the system and compute any auxiliary fields
277  Volta.Solve();
278 
279  // Determine the current size of the linear system
280  int prob_size = Volta.GetProblemSize();
281 
282  // Write fields to disk for VisIt
283  if ( visit )
284  {
285  Volta.WriteVisItFields(it);
286  }
287 
288  // Send the solution by socket to a GLVis server.
289  if (visualization)
290  {
291  Volta.DisplayToGLVis();
292  }
293  if (myid == 0 && (visit || visualization)) { cout << "done." << endl; }
294 
295  if (myid == 0)
296  {
297  cout << "AMR iteration " << it << " complete." << endl;
298  }
299 
300  // Check stopping criteria
301  if (prob_size > max_dofs)
302  {
303  if (myid == 0)
304  {
305  cout << "Reached maximum number of dofs, exiting..." << endl;
306  }
307  break;
308  }
309 
310  // Wait for user input. Ask every 10th iteration.
311  char c = 'c';
312  if (myid == 0 && (it % 10 == 0))
313  {
314  cout << "press (q)uit or (c)ontinue --> " << flush;
315  cin >> c;
316  }
317  MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
318 
319  if (c != 'c')
320  {
321  break;
322  }
323 
324  // Estimate element errors using the Zienkiewicz-Zhu error estimator.
325  Vector errors(pmesh.GetNE());
326  Volta.GetErrorEstimates(errors);
327 
328  double local_max_err = errors.Max();
329  double global_max_err;
330  MPI_Allreduce(&local_max_err, &global_max_err, 1,
331  MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
332 
333  // Make a list of elements whose error is larger than a fraction
334  // of the maximum element error. These elements will be refined.
335  Array<int> ref_list;
336  const double frac = 0.7;
337  double threshold = frac * global_max_err;
338  for (int i = 0; i < errors.Size(); i++)
339  {
340  if (errors[i] >= threshold) { ref_list.Append(i); }
341  }
342 
343  // Refine the selected elements. Since we are going to transfer the
344  // grid function x from the coarse mesh to the new fine mesh in the
345  // next step, we need to request the "two-level state" of the mesh.
346  if (myid == 0) { cout << " Refinement ..." << flush; }
347  pmesh.GeneralRefinement(ref_list);
348 
349  // Update the electrostatic solver to reflect the new state of the mesh.
350  Volta.Update();
351  }
352 
353  MPI_Finalize();
354 
355  return 0;
356 }
357 
358 // Print the Volta ascii logo to the given ostream
359 void display_banner(ostream & os)
360 {
361  os << " ____ ____ __ __ " << endl
362  << " \\ \\ / /___ | |_/ |______ " << endl
363  << " \\ Y / _ \\| |\\ __\\__ \\ " << endl
364  << " \\ ( <_> ) |_| | / __ \\_ " << endl
365  << " \\___/ \\____/|____/__| (____ / " << endl
366  << " \\/ " << endl << flush;
367 }
368 
369 // A sphere with constant permittivity. The sphere has a radius,
370 // center, and permittivity specified on the command line and stored
371 // in ds_params_.
372 double dielectric_sphere(const Vector &x)
373 {
374  double r2 = 0.0;
375 
376  for (int i=0; i<x.Size(); i++)
377  {
378  r2 += (x(i)-ds_params_(i))*(x(i)-ds_params_(i));
379  }
380 
381  if ( sqrt(r2) <= ds_params_(x.Size()) )
382  {
383  return ds_params_(x.Size()+1) * epsilon0_;
384  }
385  return epsilon0_;
386 }
387 
388 // A sphere with constant charge density. The sphere has a radius,
389 // center, and total charge specified on the command line and stored
390 // in cs_params_.
391 double charged_sphere(const Vector &x)
392 {
393  double r2 = 0.0;
394  double rho = 0.0;
395 
396  if ( cs_params_(x.Size()) > 0.0 )
397  {
398  switch ( x.Size() )
399  {
400  case 2:
401  rho = cs_params_(x.Size()+1)/(M_PI*pow(cs_params_(x.Size()),2));
402  break;
403  case 3:
404  rho = 0.75*cs_params_(x.Size()+1)/(M_PI*pow(cs_params_(x.Size()),3));
405  break;
406  default:
407  rho = 0.0;
408  }
409  }
410 
411  for (int i=0; i<x.Size(); i++)
412  {
413  r2 += (x(i)-cs_params_(i))*(x(i)-cs_params_(i));
414  }
415 
416  if ( sqrt(r2) <= cs_params_(x.Size()) )
417  {
418  return rho;
419  }
420  return 0.0;
421 }
422 
423 // A Cylindrical Rod of constant polarization. The cylinder has two
424 // axis end points, a radius, and a constant electric polarization oriented
425 // along the axis.
426 void voltaic_pile(const Vector &x, Vector &p)
427 {
428  p.SetSize(x.Size());
429  p = 0.0;
430 
431  Vector a(x.Size()); // Normalized Axis vector
432  Vector xu(x.Size()); // x vector relative to the axis end-point
433 
434  xu = x;
435 
436  for (int i=0; i<x.Size(); i++)
437  {
438  xu[i] -= vp_params_[i];
439  a[i] = vp_params_[x.Size()+i] - vp_params_[i];
440  }
441 
442  double h = a.Norml2();
443 
444  if ( h == 0.0 )
445  {
446  return;
447  }
448 
449  double r = vp_params_[2*x.Size()];
450  double xa = xu*a;
451 
452  if ( h > 0.0 )
453  {
454  xu.Add(-xa/(h*h),a);
455  }
456 
457  double xp = xu.Norml2();
458 
459  if ( xa >= 0.0 && xa <= h*h && xp <= r )
460  {
461  p.Add(vp_params_[2*x.Size()+1]/h,a);
462  }
463 }
464 
465 // To produce a uniform electric field the potential can be set
466 // to (- Ex x - Ey y - Ez z).
467 double phi_bc_uniform(const Vector &x)
468 {
469  double phi = 0.0;
470 
471  for (int i=0; i<x.Size(); i++)
472  {
473  phi -= x(i) * e_uniform_(i);
474  }
475 
476  return phi;
477 }
double phi_bc_uniform(const Vector &)
Definition: volta.cpp:467
int Size() const
Logical size of the array.
Definition: array.hpp:109
void SetSize(int s)
Resizes the vector if the new size is different.
Definition: vector.hpp:259
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:644
int Size() const
Returns the size of the vector.
Definition: vector.hpp:84
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:454
void RegisterVisItFields(VisItDataCollection &visit_dc)
void voltaic_pile(const Vector &, Vector &)
Definition: volta.cpp:426
void GetErrorEstimates(Vector &errors)
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:368
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7316
Data collection with VisIt I/O routines.
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3736
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
double charged_sphere(const Vector &)
Definition: volta.cpp:391
int SpaceDimension() const
Definition: mesh.hpp:476
int main(int argc, char *argv[])
Definition: ex1.cpp:45
MPI_Comm GetComm() const
Definition: pmesh.hpp:86
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
NURBSExtension * NURBSext
Definition: mesh.hpp:142
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:200
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
void display_banner(ostream &os)
Definition: tesla.cpp:350
void EnsureNCMesh()
Definition: mesh.cpp:6929
double dielectric_sphere(const Vector &)
Definition: volta.cpp:372
Vector data type.
Definition: vector.hpp:33
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:6864
Class for parallel meshes.
Definition: pmesh.hpp:28
bool Good() const
Definition: optparser.hpp:120