1 // Copyright (c) 2010-2021, 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
74 Coefficient * SetupPermittivityCoefficient();
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
main(int argc,char * argv[])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);
125 args.AddOption(&mesh_file, "-m", "--mesh",
126 "Mesh file to use.");
127 args.AddOption(&order, "-o", "--order",
128 "Finite element order (polynomial degree).");
129 args.AddOption(&serial_ref_levels, "-rs", "--serial-ref-levels",
130 "Number of serial refinement levels.");
131 args.AddOption(¶llel_ref_levels, "-rp", "--parallel-ref-levels",
132 "Number of parallel refinement levels.");
133 args.AddOption(&e_uniform_, "-uebc", "--uniform-e-bc",
134 "Specify if the three components of the constant "
135 "electric field");
136 args.AddOption(&pw_eps_, "-pwe", "--piecewise-eps",
137 "Piecewise values of Permittivity");
138 args.AddOption(&ds_params_, "-ds", "--dielectric-sphere-params",
139 "Center, Radius, and Permittivity of Dielectric Sphere");
140 args.AddOption(&cs_params_, "-cs", "--charged-sphere-params",
141 "Center, Radius, and Total Charge of Charged Sphere");
142 args.AddOption(&pc_params_, "-pc", "--point-charge-params",
143 "Charges and locations of Point Charges");
144 args.AddOption(&vp_params_, "-vp", "--voltaic-pile-params",
145 "Axis End Points, Radius, and "
146 "Polarization of Cylindrical Voltaic Pile");
147 args.AddOption(&dbcs, "-dbcs", "--dirichlet-bc-surf",
148 "Dirichlet Boundary Condition Surfaces");
149 args.AddOption(&dbcv, "-dbcv", "--dirichlet-bc-vals",
150 "Dirichlet Boundary Condition Values");
151 args.AddOption(&dbcg, "-dbcg", "--dirichlet-bc-gradient",
152 "-no-dbcg", "--no-dirichlet-bc-gradient",
153 "Dirichlet Boundary Condition Gradient (phi = -z)");
154 args.AddOption(&nbcs, "-nbcs", "--neumann-bc-surf",
155 "Neumann Boundary Condition Surfaces");
156 args.AddOption(&nbcv, "-nbcv", "--neumann-bc-vals",
157 "Neumann Boundary Condition Values");
158 args.AddOption(&maxit, "-maxit", "--max-amr-iterations",
159 "Max number of iterations in the main AMR loop.");
160 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
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
249 Coefficient * epsCoef = SetupPermittivityCoefficient();
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
display_banner(ostream & os)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 *
SetupPermittivityCoefficient()391 SetupPermittivityCoefficient()
392 {
393 Coefficient * coef = NULL;
394
395 if ( ds_params_.Size() > 0 )
396 {
397 coef = new FunctionCoefficient(dielectric_sphere);
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_.
dielectric_sphere(const Vector & x)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_.
charged_sphere(const Vector & x)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.
voltaic_pile(const Vector & x,Vector & p)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).
phi_bc_uniform(const Vector & x)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 }
522