MFEM  v3.3
Finite element discretization library
tesla.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 // Tesla Miniapp: Simple Magnetostatics Simulation Code
14 // -----------------------------------------------------
15 //
16 // This miniapp solves a simple 3D magnetostatic problem.
17 //
18 // Curl 1/mu Curl A = J + Curl mu0/mu M
19 //
20 // The permeability function is that of the vacuum with an optional diamagnetic
21 // or paramagnetic spherical shell. The optional current density takes the form
22 // of a user defined ring of current. The optional magnetization consists of a
23 // cylindrical bar of constant magnetization.
24 //
25 // The boundary conditions either apply a user selected uniform magnetic flux
26 // density or a surface current flowing between user defined surfaces.
27 //
28 // We discretize the vector potential with H(Curl) finite elements. The magnetic
29 // flux B is discretized with H(Div) finite elements.
30 //
31 // Compile with: make tesla
32 //
33 // Sample runs:
34 //
35 // A cylindrical bar magnet in a metal sphere:
36 // mpirun -np 4 tesla -bm '0 -0.5 0 0 0.5 0 0.2 1'
37 //
38 // A spherical shell of paramagnetic material in a uniform B field:
39 // mpirun -np 4 tesla -ubbc '0 0 1' -ms '0 0 0 0.2 0.4 10'
40 //
41 // A ring of current in a metal sphere:
42 // mpirun -np 4 tesla -cr '0 0 -0.2 0 0 0.2 0.2 0.4 1'
43 //
44 // A Halbach array of permanent magnets:
45 // mpirun -np 4 tesla -m ../../data/beam-hex.mesh -rs 2
46 // -ha '1 0.1 0.3 7 0.9 0.7 0 1 12'
47 //
48 // An example demonstrating the use of surface currents:
49 // mpirun -np 4 tesla -m square-angled-pipe.mesh
50 // -kbcs '3' -vbcs '1 2' -vbcv '-0.5 0.5'
51 //
52 // An example combining the paramagnetic shell, permanent magnet,
53 // and current ring:
54 // mpirun -np 4 tesla -m ../../data/inline-hex.mesh
55 // -ms '0.5 0.5 0.5 0.4 0.45 20'
56 // -bm '0.5 0.5 0.3 0.5 0.5 0.7 0.1 1'
57 // -cr '0.5 0.5 0.45 0.5 0.5 0.55 0.2 0.3 1'
58 //
59 // By default the sources and fields are all zero:
60 // mpirun -np 4 tesla
61 
62 #include "tesla_solver.hpp"
63 #include <fstream>
64 #include <iostream>
65 
66 using namespace std;
67 using namespace mfem;
68 using namespace mfem::electromagnetics;
69 
70 // Permeability Function
72 
73 static Vector pw_mu_(0); // Piecewise permeability values
74 static Vector pw_mu_inv_(0); // Piecewise inverse permeability values
75 static Vector ms_params_(0); // Center, Inner and Outer Radii, and
76 // Permeability of magnetic shell
77 double magnetic_shell(const Vector &);
78 double magnetic_shell_inv(const Vector & x) { return 1.0/magnetic_shell(x); }
79 
80 // Current Density Function
81 static Vector cr_params_(0); // Axis Start, Axis End, Inner Ring Radius,
82 // Outer Ring Radius, and Total Current
83 // of current ring (annulus)
84 void current_ring(const Vector &, Vector &);
85 
86 // Magnetization
87 static Vector bm_params_(0); // Axis Start, Axis End, Bar Radius,
88 // and Magnetic Field Magnitude
89 void bar_magnet(const Vector &, Vector &);
90 
91 static Vector ha_params_(0); // Bounding box,
92 // axis index (0->'x', 1->'y', 2->'z'),
93 // rotation axis index
94 // and number of segments
95 void halbach_array(const Vector &, Vector &);
96 
97 // A Field Boundary Condition for B = (Bx,By,Bz)
98 static Vector b_uniform_(0);
99 void a_bc_uniform(const Vector &, Vector&);
100 
101 // Phi_M Boundary Condition for H = (0,0,1)
102 double phi_m_bc_uniform(const Vector &x);
103 
104 // Prints the program's logo to the given output stream
105 void display_banner(ostream & os);
106 
107 int main(int argc, char *argv[])
108 {
109  MPI_Session mpi(argc, argv);
110 
111  if ( mpi.Root() ) { display_banner(cout); }
112 
113  // Parse command-line options.
114  const char *mesh_file = "../../data/ball-nurbs.mesh";
115  int order = 1;
116  int maxit = 100;
117  int serial_ref_levels = 0;
118  int parallel_ref_levels = 0;
119  bool visualization = true;
120  bool visit = true;
121 
122  Array<int> kbcs;
123  Array<int> vbcs;
124 
125  Vector vbcv;
126 
127  OptionsParser args(argc, argv);
128  args.AddOption(&mesh_file, "-m", "--mesh",
129  "Mesh file to use.");
130  args.AddOption(&order, "-o", "--order",
131  "Finite element order (polynomial degree).");
132  args.AddOption(&serial_ref_levels, "-rs", "--serial-ref-levels",
133  "Number of serial refinement levels.");
134  args.AddOption(&parallel_ref_levels, "-rp", "--parallel-ref-levels",
135  "Number of parallel refinement levels.");
136  args.AddOption(&b_uniform_, "-ubbc", "--uniform-b-bc",
137  "Specify if the three components of the constant magnetic flux density");
138  args.AddOption(&pw_mu_, "-pwm", "--piecewise-mu",
139  "Piecewise values of Permeability");
140  args.AddOption(&ms_params_, "-ms", "--magnetic-shell-params",
141  "Center, Inner Radius, Outer Radius, and Permeability of Magnetic Shell");
142  args.AddOption(&cr_params_, "-cr", "--current-ring-params",
143  "Axis End Points, Inner Radius, Outer Radius and Total Current of Annulus");
144  args.AddOption(&bm_params_, "-bm", "--bar-magnet-params",
145  "Axis End Points, Radius, and Magnetic Field of Cylindrical Magnet");
146  args.AddOption(&ha_params_, "-ha", "--halbach-array-params",
147  "Bounding Box Corners and Number of Segments");
148  args.AddOption(&kbcs, "-kbcs", "--surface-current-bc",
149  "Surfaces for the Surface Current (K) Boundary Condition");
150  args.AddOption(&vbcs, "-vbcs", "--voltage-bc-surf",
151  "Voltage Boundary Condition Surfaces (to drive K)");
152  args.AddOption(&vbcv, "-vbcv", "--voltage-bc-vals",
153  "Voltage Boundary Condition Values (to drive K)");
154  args.AddOption(&maxit, "-maxit", "--max-amr-iterations",
155  "Max number of iterations in the main AMR loop.");
156  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
157  "--no-visualization",
158  "Enable or disable GLVis visualization.");
159  args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit",
160  "Enable or disable VisIt visualization.");
161  args.Parse();
162  if (!args.Good())
163  {
164  if (mpi.Root())
165  {
166  args.PrintUsage(cout);
167  }
168  return 1;
169  }
170  if (mpi.Root())
171  {
172  args.PrintOptions(cout);
173  }
174 
175  // Read the (serial) mesh from the given mesh file on all processors. We
176  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
177  // and volume meshes with the same code.
178  Mesh *mesh = new Mesh(mesh_file, 1, 1);
179 
180  if (mpi.Root())
181  {
182  cout << "Starting initialization." << endl;
183  }
184 
185  // Project a NURBS mesh to a piecewise-quadratic curved mesh
186  if (mesh->NURBSext)
187  {
188  mesh->UniformRefinement();
189  if (serial_ref_levels > 0) { serial_ref_levels--; }
190 
191  mesh->SetCurvature(2);
192  }
193 
194  // Ensure that quad and hex meshes are treated as non-conforming.
195  mesh->EnsureNCMesh();
196 
197  // Refine the serial mesh on all processors to increase the resolution. In
198  // this example we do 'ref_levels' of uniform refinement.
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 values for Voltage BCs were not set issue a warning and exit
218  if ( ( vbcs.Size() > 0 && kbcs.Size() == 0 ) ||
219  ( kbcs.Size() > 0 && vbcs.Size() == 0 ) ||
220  ( vbcv.Size() < vbcs.Size() ) )
221  {
222  if ( mpi.Root() )
223  {
224  cout << "The surface current (K) boundary condition requires "
225  << "surface current boundary condition surfaces (with -kbcs), "
226  << "voltage boundary condition surface (with -vbcs), "
227  << "and voltage boundary condition values (with -vbcv)."
228  << endl;
229  }
230  return 3;
231  }
232 
233  // Create a coefficient describing the magnetic permeability
235 
236  // Create the Magnetostatic solver
237  TeslaSolver Tesla(pmesh, order, kbcs, vbcs, vbcv, *muInvCoef,
238  (b_uniform_.Size() > 0 ) ? a_bc_uniform : NULL,
239  (cr_params_.Size() > 0 ) ? current_ring : NULL,
240  (bm_params_.Size() > 0 ) ? bar_magnet :
241  (ha_params_.Size() > 0 ) ? halbach_array : NULL);
242 
243  // Initialize GLVis visualization
244  if (visualization)
245  {
246  Tesla.InitializeGLVis();
247  }
248 
249  // Initialize VisIt visualization
250  VisItDataCollection visit_dc("Tesla-AMR-Parallel", &pmesh);
251 
252  if ( visit )
253  {
254  Tesla.RegisterVisItFields(visit_dc);
255  }
256  if (mpi.Root()) { cout << "Initialization done." << endl; }
257 
258  // The main AMR loop. In each iteration we solve the problem on the current
259  // mesh, visualize the solution, estimate the error on all elements, refine
260  // the worst elements and update all objects to work with the new mesh. We
261  // refine until the maximum number of dofs in the Nedelec finite element
262  // space reaches 10 million.
263  const int max_dofs = 10000000;
264  for (int it = 1; it <= maxit; it++)
265  {
266  if (mpi.Root())
267  {
268  cout << "\nAMR Iteration " << it << endl;
269  }
270 
271  // Display the current number of DoFs in each finite element space
272  Tesla.PrintSizes();
273 
274  // Assemble all forms
275  Tesla.Assemble();
276 
277  // Solve the system and compute any auxiliary fields
278  Tesla.Solve();
279 
280  // Determine the current size of the linear system
281  int prob_size = Tesla.GetProblemSize();
282 
283  // Write fields to disk for VisIt
284  if ( visit )
285  {
286  Tesla.WriteVisItFields(it);
287  }
288 
289  // Send the solution by socket to a GLVis server.
290  if (visualization)
291  {
292  Tesla.DisplayToGLVis();
293  }
294 
295  if (mpi.Root())
296  {
297  cout << "AMR iteration " << it << " complete." << endl;
298  }
299 
300  // Check stopping criteria
301  if (prob_size > max_dofs)
302  {
303  if (mpi.Root())
304  {
305  cout << "Reached maximum number of dofs, exiting..." << endl;
306  }
307  break;
308  }
309  if ( it == maxit )
310  {
311  break;
312  }
313 
314  // Wait for user input. Ask every 10th iteration.
315  char c = 'c';
316  if (mpi.Root() && (it % 10 == 0))
317  {
318  cout << "press (q)uit or (c)ontinue --> " << flush;
319  cin >> c;
320  }
321  MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
322 
323  if (c != 'c')
324  {
325  break;
326  }
327 
328  // Estimate element errors using the Zienkiewicz-Zhu error estimator.
329  Vector errors(pmesh.GetNE());
330  Tesla.GetErrorEstimates(errors);
331 
332  double local_max_err = errors.Max();
333  double global_max_err;
334  MPI_Allreduce(&local_max_err, &global_max_err, 1,
335  MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
336 
337  // Refine the elements whose error is larger than a fraction of the
338  // maximum element error.
339  const double frac = 0.5;
340  double threshold = frac * global_max_err;
341  if (mpi.Root()) { cout << "Refining ..." << endl; }
342  pmesh.RefineByError(errors, threshold);
343 
344  // Update the magnetostatic solver to reflect the new state of the mesh.
345  Tesla.Update();
346 
347  if (pmesh.Nonconforming() && mpi.WorldSize() > 1)
348  {
349  if (mpi.Root()) { cout << "Rebalancing ..." << endl; }
350  pmesh.Rebalance();
351 
352  // Update again after rebalancing
353  Tesla.Update();
354  }
355  }
356 
357  delete muInvCoef;
358 
359  return 0;
360 }
361 
362 // Print the Volta ascii logo to the given ostream
363 void display_banner(ostream & os)
364 {
365  os << " ___________ __ " << endl
366  << " \\__ ___/___ _____| | _____ " << endl
367  << " | |_/ __ \\ / ___/ | \\__ \\ " << endl
368  << " | |\\ ___/ \\___ \\| |__/ __ \\_ " << endl
369  << " |____| \\___ >____ >____(____ / " << endl
370  << " \\/ \\/ \\/ " << endl << flush;
371 }
372 
373 // The Permeability is a required coefficient which may be defined in
374 // various ways so we'll determine the appropriate coefficient type here.
375 Coefficient *
377 {
378  Coefficient * coef = NULL;
379 
380  if ( ms_params_.Size() > 0 )
381  {
383  }
384  else if ( pw_mu_.Size() > 0 )
385  {
386  pw_mu_inv_.SetSize(pw_mu_.Size());
387  for (int i = 0; i < pw_mu_.Size(); i++)
388  {
389  MFEM_ASSERT( pw_mu_[i] > 0.0, "permeability values must be positive" );
390  pw_mu_inv_[i] = 1.0/pw_mu_[i];
391  }
392  coef = new PWConstCoefficient(pw_mu_inv_);
393  }
394  else
395  {
396  coef = new ConstantCoefficient(1.0/mu0_);
397  }
398 
399  return coef;
400 }
401 
402 // A spherical shell with constant permeability. The sphere has inner
403 // and outer radii, center, and relative permeability specified on the
404 // command line and stored in ms_params_.
405 double magnetic_shell(const Vector &x)
406 {
407  double r2 = 0.0;
408 
409  for (int i = 0; i < x.Size(); i++)
410  {
411  r2 += (x(i) - ms_params_(i))*(x(i) - ms_params_(i));
412  }
413 
414  if ( sqrt(r2) >= ms_params_(x.Size()) &&
415  sqrt(r2) <= ms_params_(x.Size()+1) )
416  {
417  return mu0_*ms_params_(x.Size()+2);
418  }
419  return mu0_;
420 }
421 
422 // An annular ring of current density. The ring has two axis end
423 // points, inner and outer radii, and a constant current in Amperes.
424 void current_ring(const Vector &x, Vector &j)
425 {
426  MFEM_ASSERT(x.Size() == 3, "current_ring source requires 3D space.");
427 
428  j.SetSize(x.Size());
429  j = 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  Vector ju(x.Size()); // Unit vector in direction of current
434 
435  xu = x;
436 
437  for (int i=0; i<x.Size(); i++)
438  {
439  xu[i] -= cr_params_[i];
440  a[i] = cr_params_[x.Size()+i] - cr_params_[i];
441  }
442 
443  double h = a.Norml2();
444 
445  if ( h == 0.0 )
446  {
447  return;
448  }
449 
450  double ra = cr_params_[2*x.Size()+0];
451  double rb = cr_params_[2*x.Size()+1];
452  if ( ra > rb )
453  {
454  double rc = ra;
455  ra = rb;
456  rb = rc;
457  }
458  double xa = xu*a;
459 
460  if ( h > 0.0 )
461  {
462  xu.Add(-xa/(h*h),a);
463  }
464 
465  double xp = xu.Norml2();
466 
467  if ( xa >= 0.0 && xa <= h*h && xp >= ra && xp <= rb )
468  {
469  ju(0) = a(1) * xu(2) - a(2) * xu(1);
470  ju(1) = a(2) * xu(0) - a(0) * xu(2);
471  ju(2) = a(0) * xu(1) - a(1) * xu(0);
472  ju /= h;
473 
474  j.Add(cr_params_[2*x.Size()+2]/(h*(rb-ra)),ju);
475  }
476 }
477 
478 // A Cylindrical Rod of constant magnetization. The cylinder has two
479 // axis end points, a radius, and a constant magnetic field oriented
480 // along the axis.
481 void bar_magnet(const Vector &x, Vector &m)
482 {
483  m.SetSize(x.Size());
484  m = 0.0;
485 
486  Vector a(x.Size()); // Normalized Axis vector
487  Vector xu(x.Size()); // x vector relative to the axis end-point
488 
489  xu = x;
490 
491  for (int i=0; i<x.Size(); i++)
492  {
493  xu[i] -= bm_params_[i];
494  a[i] = bm_params_[x.Size()+i] - bm_params_[i];
495  }
496 
497  double h = a.Norml2();
498 
499  if ( h == 0.0 )
500  {
501  return;
502  }
503 
504  double r = bm_params_[2*x.Size()];
505  double xa = xu*a;
506 
507  if ( h > 0.0 )
508  {
509  xu.Add(-xa/(h*h),a);
510  }
511 
512  double xp = xu.Norml2();
513 
514  if ( xa >= 0.0 && xa <= h*h && xp <= r )
515  {
516  m.Add(bm_params_[2*x.Size()+1]/h,a);
517  }
518 }
519 
520 // A Square Rod of rotating magnetized segments. The rod is defined
521 // by a bounding box and a number of segments. The magnetization in
522 // each segment is constant and follows a rotating pattern.
523 void halbach_array(const Vector &x, Vector &m)
524 {
525  m.SetSize(x.Size());
526  m = 0.0;
527 
528  // Check Bounding Box
529  if ( x[0] < ha_params_[0] || x[0] > ha_params_[3] ||
530  x[1] < ha_params_[1] || x[1] > ha_params_[4] ||
531  x[2] < ha_params_[2] || x[2] > ha_params_[5] )
532  {
533  return;
534  }
535 
536  int ai = (int)ha_params_[6];
537  int ri = (int)ha_params_[7];
538  int n = (int)ha_params_[8];
539 
540  int i = (int)n * (x[ai] - ha_params_[ai]) /
541  (ha_params_[ai+3] - ha_params_[ai]);
542 
543  m[(ri + 1 + (i % 2)) % 3] = pow(-1.0,i/2);
544 }
545 
546 // To produce a uniform magnetic flux the vector potential can be set
547 // to ( By z, Bz x, Bx y).
548 void a_bc_uniform(const Vector & x, Vector & a)
549 {
550  a.SetSize(3);
551  a(0) = b_uniform_(1) * x(2);
552  a(1) = b_uniform_(2) * x(0);
553  a(2) = b_uniform_(0) * x(1);
554 }
555 
556 // To produce a uniform magnetic field the scalar potential can be set
557 // to -z (or -y in 2D).
558 double phi_m_bc_uniform(const Vector &x)
559 {
560  return -x(x.Size()-1);
561 }
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
double magnetic_shell(const Vector &)
Definition: tesla.cpp:405
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
void bar_magnet(const Vector &, Vector &)
Definition: tesla.cpp:481
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:584
void GetErrorEstimates(Vector &errors)
STL namespace.
double phi_m_bc_uniform(const Vector &x)
Definition: tesla.cpp:558
void Rebalance()
Load balance the mesh. NC meshes only.
Definition: pmesh.cpp:2713
bool Nonconforming() const
Definition: mesh.hpp:968
A simple convenience class that calls MPI_Init() at construction and MPI_Finalize() at destruction...
double magnetic_shell_inv(const Vector &x)
Definition: tesla.cpp:78
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6684
void current_ring(const Vector &, Vector &)
Definition: tesla.cpp:424
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.
int main(int argc, char *argv[])
Definition: tesla.cpp:107
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void RegisterVisItFields(VisItDataCollection &visit_dc)
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
void a_bc_uniform(const Vector &, Vector &)
Definition: tesla.cpp:548
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
Coefficient * SetupInvPermeabilityCoefficient()
Definition: tesla.cpp:376
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
void display_banner(ostream &os)
Definition: tesla.cpp:363
void EnsureNCMesh(bool triangles_nonconforming=false)
Definition: mesh.cpp:6350
class for C-function coefficient
Vector data type.
Definition: vector.hpp:36
Class for parallel meshes.
Definition: pmesh.hpp:28
bool Good() const
Definition: optparser.hpp:120
void halbach_array(const Vector &, Vector &)
Definition: tesla.cpp:523