MFEM  v4.4.0 Finite element discretization library
tesla.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 // 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::Init(argc, argv);
110  Hypre::Init();
111
112  if ( Mpi::Root() ) { display_banner(cout); }
113
114  // Parse command-line options.
115  const char *mesh_file = "../../data/ball-nurbs.mesh";
116  int order = 1;
117  int maxit = 100;
118  int serial_ref_levels = 0;
119  int parallel_ref_levels = 0;
120  bool visualization = true;
121  bool visit = true;
122
123  Array<int> kbcs;
124  Array<int> vbcs;
125
126  Vector vbcv;
127
128  OptionsParser args(argc, argv);
130  "Mesh file to use.");
132  "Finite element order (polynomial degree).");
134  "Number of serial refinement levels.");
136  "Number of parallel refinement levels.");
138  "Specify if the three components of the constant magnetic flux density");
140  "Piecewise values of Permeability");
144  "Axis End Points, Inner Radius, Outer Radius and Total Current of Annulus");
146  "Axis End Points, Radius, and Magnetic Field of Cylindrical Magnet");
148  "Bounding Box Corners and Number of Segments");
150  "Surfaces for the Surface Current (K) Boundary Condition");
152  "Voltage Boundary Condition Surfaces (to drive K)");
154  "Voltage Boundary Condition Values (to drive K)");
156  "Max number of iterations in the main AMR loop.");
158  "--no-visualization",
159  "Enable or disable GLVis visualization.");
160  args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit",
161  "Enable or disable VisIt visualization.");
162  args.Parse();
163  if (!args.Good())
164  {
165  if (Mpi::Root())
166  {
167  args.PrintUsage(cout);
168  }
169  return 1;
170  }
171  if (Mpi::Root())
172  {
173  args.PrintOptions(cout);
174  }
175
176  // Read the (serial) mesh from the given mesh file on all processors. We
177  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
178  // and volume meshes with the same code.
179  Mesh *mesh = new Mesh(mesh_file, 1, 1);
180
181  if (Mpi::Root())
182  {
183  cout << "Starting initialization." << endl;
184  }
185
186  // Project a NURBS mesh to a piecewise-quadratic curved mesh
187  if (mesh->NURBSext)
188  {
189  mesh->UniformRefinement();
190  if (serial_ref_levels > 0) { serial_ref_levels--; }
191
192  mesh->SetCurvature(2);
193  }
194
195  // Ensure that quad and hex meshes are treated as non-conforming.
196  mesh->EnsureNCMesh();
197
198  // Refine the serial mesh on all processors to increase the resolution. In
199  // this example we do 'ref_levels' of uniform refinement.
200  for (int l = 0; l < serial_ref_levels; l++)
201  {
202  mesh->UniformRefinement();
203  }
204
205  // Define a parallel mesh by a partitioning of the serial mesh. Refine
206  // this mesh further in parallel to increase the resolution. Once the
207  // parallel mesh is defined, the serial mesh can be deleted.
208  ParMesh pmesh(MPI_COMM_WORLD, *mesh);
209  delete mesh;
210
211  // Refine this mesh in parallel to increase the resolution.
212  int par_ref_levels = parallel_ref_levels;
213  for (int l = 0; l < par_ref_levels; l++)
214  {
215  pmesh.UniformRefinement();
216  }
217  // Make sure tet-only meshes are marked for local refinement.
218  pmesh.Finalize(true);
219
220  // If values for Voltage BCs were not set issue a warning and exit
221  if ( ( vbcs.Size() > 0 && kbcs.Size() == 0 ) ||
222  ( kbcs.Size() > 0 && vbcs.Size() == 0 ) ||
223  ( vbcv.Size() < vbcs.Size() ) )
224  {
225  if ( Mpi::Root() )
226  {
227  cout << "The surface current (K) boundary condition requires "
228  << "surface current boundary condition surfaces (with -kbcs), "
229  << "voltage boundary condition surface (with -vbcs), "
230  << "and voltage boundary condition values (with -vbcv)."
231  << endl;
232  }
233  return 3;
234  }
235
236  // Create a coefficient describing the magnetic permeability
238
239  // Create the Magnetostatic solver
240  TeslaSolver Tesla(pmesh, order, kbcs, vbcs, vbcv, *muInvCoef,
241  (b_uniform_.Size() > 0 ) ? a_bc_uniform : NULL,
242  (cr_params_.Size() > 0 ) ? current_ring : NULL,
243  (bm_params_.Size() > 0 ) ? bar_magnet :
244  (ha_params_.Size() > 0 ) ? halbach_array : NULL);
245
246  // Initialize GLVis visualization
247  if (visualization)
248  {
249  Tesla.InitializeGLVis();
250  }
251
252  // Initialize VisIt visualization
253  VisItDataCollection visit_dc("Tesla-AMR-Parallel", &pmesh);
254
255  if ( visit )
256  {
257  Tesla.RegisterVisItFields(visit_dc);
258  }
259  if (Mpi::Root()) { cout << "Initialization done." << endl; }
260
261  // The main AMR loop. In each iteration we solve the problem on the current
262  // mesh, visualize the solution, estimate the error on all elements, refine
263  // the worst elements and update all objects to work with the new mesh. We
264  // refine until the maximum number of dofs in the Nedelec finite element
265  // space reaches 10 million.
266  const int max_dofs = 10000000;
267  for (int it = 1; it <= maxit; it++)
268  {
269  if (Mpi::Root())
270  {
271  cout << "\nAMR Iteration " << it << endl;
272  }
273
274  // Display the current number of DoFs in each finite element space
275  Tesla.PrintSizes();
276
277  // Assemble all forms
278  Tesla.Assemble();
279
280  // Solve the system and compute any auxiliary fields
281  Tesla.Solve();
282
283  // Determine the current size of the linear system
284  int prob_size = Tesla.GetProblemSize();
285
286  // Write fields to disk for VisIt
287  if ( visit )
288  {
289  Tesla.WriteVisItFields(it);
290  }
291
292  // Send the solution by socket to a GLVis server.
293  if (visualization)
294  {
295  Tesla.DisplayToGLVis();
296  }
297
298  if (Mpi::Root())
299  {
300  cout << "AMR iteration " << it << " complete." << endl;
301  }
302
303  // Check stopping criteria
304  if (prob_size > max_dofs)
305  {
306  if (Mpi::Root())
307  {
308  cout << "Reached maximum number of dofs, exiting..." << endl;
309  }
310  break;
311  }
312  if ( it == maxit )
313  {
314  break;
315  }
316
317  // Wait for user input. Ask every 10th iteration.
318  char c = 'c';
319  if (Mpi::Root() && (it % 10 == 0))
320  {
321  cout << "press (q)uit or (c)ontinue --> " << flush;
322  cin >> c;
323  }
324  MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
325
326  if (c != 'c')
327  {
328  break;
329  }
330
331  // Estimate element errors using the Zienkiewicz-Zhu error estimator.
332  Vector errors(pmesh.GetNE());
333  Tesla.GetErrorEstimates(errors);
334
335  double local_max_err = errors.Max();
336  double global_max_err;
337  MPI_Allreduce(&local_max_err, &global_max_err, 1,
338  MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
339
340  // Refine the elements whose error is larger than a fraction of the
341  // maximum element error.
342  const double frac = 0.5;
343  double threshold = frac * global_max_err;
344  if (Mpi::Root()) { cout << "Refining ..." << endl; }
345  pmesh.RefineByError(errors, threshold);
346
347  // Update the magnetostatic solver to reflect the new state of the mesh.
348  Tesla.Update();
349
350  if (pmesh.Nonconforming() && Mpi::WorldSize() > 1)
351  {
352  if (Mpi::Root()) { cout << "Rebalancing ..." << endl; }
353  pmesh.Rebalance();
354
355  // Update again after rebalancing
356  Tesla.Update();
357  }
358  }
359
360  delete muInvCoef;
361
362  return 0;
363 }
364
365 // Print the Volta ascii logo to the given ostream
366 void display_banner(ostream & os)
367 {
368  os << " ___________ __ " << endl
369  << " \\__ ___/___ _____| | _____ " << endl
370  << " | |_/ __ \\ / ___/ | \\__ \\ " << endl
371  << " | |\\ ___/ \\___ \\| |__/ __ \\_ " << endl
372  << " |____| \\___ >____ >____(____ / " << endl
373  << " \\/ \\/ \\/ " << endl << flush;
374 }
375
376 // The Permeability is a required coefficient which may be defined in
377 // various ways so we'll determine the appropriate coefficient type here.
378 Coefficient *
380 {
381  Coefficient * coef = NULL;
382
383  if ( ms_params_.Size() > 0 )
384  {
386  }
387  else if ( pw_mu_.Size() > 0 )
388  {
389  pw_mu_inv_.SetSize(pw_mu_.Size());
390  for (int i = 0; i < pw_mu_.Size(); i++)
391  {
392  MFEM_ASSERT( pw_mu_[i] > 0.0, "permeability values must be positive" );
393  pw_mu_inv_[i] = 1.0/pw_mu_[i];
394  }
395  coef = new PWConstCoefficient(pw_mu_inv_);
396  }
397  else
398  {
399  coef = new ConstantCoefficient(1.0/mu0_);
400  }
401
402  return coef;
403 }
404
405 // A spherical shell with constant permeability. The sphere has inner
406 // and outer radii, center, and relative permeability specified on the
407 // command line and stored in ms_params_.
408 double magnetic_shell(const Vector &x)
409 {
410  double r2 = 0.0;
411
412  for (int i = 0; i < x.Size(); i++)
413  {
414  r2 += (x(i) - ms_params_(i))*(x(i) - ms_params_(i));
415  }
416
417  if ( sqrt(r2) >= ms_params_(x.Size()) &&
418  sqrt(r2) <= ms_params_(x.Size()+1) )
419  {
420  return mu0_*ms_params_(x.Size()+2);
421  }
422  return mu0_;
423 }
424
425 // An annular ring of current density. The ring has two axis end
426 // points, inner and outer radii, and a constant current in Amperes.
427 void current_ring(const Vector &x, Vector &j)
428 {
429  MFEM_ASSERT(x.Size() == 3, "current_ring source requires 3D space.");
430
431  j.SetSize(x.Size());
432  j = 0.0;
433
434  Vector a(x.Size()); // Normalized Axis vector
435  Vector xu(x.Size()); // x vector relative to the axis end-point
436  Vector ju(x.Size()); // Unit vector in direction of current
437
438  xu = x;
439
440  for (int i=0; i<x.Size(); i++)
441  {
442  xu[i] -= cr_params_[i];
443  a[i] = cr_params_[x.Size()+i] - cr_params_[i];
444  }
445
446  double h = a.Norml2();
447
448  if ( h == 0.0 )
449  {
450  return;
451  }
452
453  double ra = cr_params_[2*x.Size()+0];
454  double rb = cr_params_[2*x.Size()+1];
455  if ( ra > rb )
456  {
457  double rc = ra;
458  ra = rb;
459  rb = rc;
460  }
461  double xa = xu*a;
462
463  if ( h > 0.0 )
464  {
466  }
467
468  double xp = xu.Norml2();
469
470  if ( xa >= 0.0 && xa <= h*h && xp >= ra && xp <= rb )
471  {
472  ju(0) = a(1) * xu(2) - a(2) * xu(1);
473  ju(1) = a(2) * xu(0) - a(0) * xu(2);
474  ju(2) = a(0) * xu(1) - a(1) * xu(0);
475  ju /= h;
476
478  }
479 }
480
481 // A Cylindrical Rod of constant magnetization. The cylinder has two
482 // axis end points, a radius, and a constant magnetic field oriented
483 // along the axis.
484 void bar_magnet(const Vector &x, Vector &m)
485 {
486  m.SetSize(x.Size());
487  m = 0.0;
488
489  Vector a(x.Size()); // Normalized Axis vector
490  Vector xu(x.Size()); // x vector relative to the axis end-point
491
492  xu = x;
493
494  for (int i=0; i<x.Size(); i++)
495  {
496  xu[i] -= bm_params_[i];
497  a[i] = bm_params_[x.Size()+i] - bm_params_[i];
498  }
499
500  double h = a.Norml2();
501
502  if ( h == 0.0 )
503  {
504  return;
505  }
506
507  double r = bm_params_[2*x.Size()];
508  double xa = xu*a;
509
510  if ( h > 0.0 )
511  {
513  }
514
515  double xp = xu.Norml2();
516
517  if ( xa >= 0.0 && xa <= h*h && xp <= r )
518  {
520  }
521 }
522
523 // A Square Rod of rotating magnetized segments. The rod is defined
524 // by a bounding box and a number of segments. The magnetization in
525 // each segment is constant and follows a rotating pattern.
526 void halbach_array(const Vector &x, Vector &m)
527 {
528  m.SetSize(x.Size());
529  m = 0.0;
530
531  // Check Bounding Box
532  if ( x[0] < ha_params_[0] || x[0] > ha_params_[3] ||
533  x[1] < ha_params_[1] || x[1] > ha_params_[4] ||
534  x[2] < ha_params_[2] || x[2] > ha_params_[5] )
535  {
536  return;
537  }
538
539  int ai = (int)ha_params_[6];
540  int ri = (int)ha_params_[7];
541  int n = (int)ha_params_[8];
542
543  int i = (int)n * (x[ai] - ha_params_[ai]) /
544  (ha_params_[ai+3] - ha_params_[ai]);
545
546  m[(ri + 1 + (i % 2)) % 3] = pow(-1.0,i/2);
547 }
548
549 // To produce a uniform magnetic flux the vector potential can be set
550 // to ( By z, Bz x, Bx y).
551 void a_bc_uniform(const Vector & x, Vector & a)
552 {
553  a.SetSize(3);
554  a(0) = b_uniform_(1) * x(2);
555  a(1) = b_uniform_(2) * x(0);
556  a(2) = b_uniform_(0) * x(1);
557 }
558
559 // To produce a uniform magnetic field the scalar potential can be set
560 // to -z (or -y in 2D).
561 double phi_m_bc_uniform(const Vector &x)
562 {
563  return -x(x.Size()-1);
564 }
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
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:521
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
void bar_magnet(const Vector &, Vector &)
Definition: tesla.cpp:484
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:928
void GetErrorEstimates(Vector &errors)
double phi_m_bc_uniform(const Vector &x)
Definition: tesla.cpp:561
void Rebalance()
Definition: pmesh.cpp:3903
bool Nonconforming() const
Definition: mesh.hpp:1626
void Finalize(bool refine=false, bool fix_orientation=false) override
Finalize the construction of a general Mesh.
Definition: pmesh.cpp:1515
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
double magnetic_shell_inv(const Vector &x)
Definition: tesla.cpp:78
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition: mesh.cpp:9394
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9737
void current_ring(const Vector &, Vector &)
Definition: tesla.cpp:427
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:5312
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void RegisterVisItFields(VisItDataCollection &visit_dc)
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
MPI_Comm GetComm() const
Definition: pmesh.hpp:288
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:9463
void a_bc_uniform(const Vector &, Vector &)
Definition: tesla.cpp:551
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:272
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:252
Coefficient * SetupInvPermeabilityCoefficient()
Definition: tesla.cpp:379
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
double magnetic_shell(const Vector &)
Definition: maxwell.cpp:401
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
void display_banner(ostream &os)
Definition: joule.cpp:778
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
void halbach_array(const Vector &, Vector &)
Definition: tesla.cpp:526