MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
tesla.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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// 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
66using namespace std;
67using namespace mfem;
68using namespace mfem::electromagnetics;
69
70// Permeability Function
72
73static Vector pw_mu_(0); // Piecewise permeability values
74static Vector pw_mu_inv_(0); // Piecewise inverse permeability values
75static Vector ms_params_(0); // Center, Inner and Outer Radii, and
76// Permeability of magnetic shell
78real_t magnetic_shell_inv(const Vector & x) { return 1.0/magnetic_shell(x); }
79
80// Current Density Function
81static Vector cr_params_(0); // Axis Start, Axis End, Inner Ring Radius,
82// Outer Ring Radius, and Total Current
83// of current ring (annulus)
84void current_ring(const Vector &, Vector &);
85
86// Magnetization
87static Vector bm_params_(0); // Axis Start, Axis End, Bar Radius,
88// and Magnetic Field Magnitude
89void bar_magnet(const Vector &, Vector &);
90
91static Vector ha_params_(0); // Bounding box,
92// axis index (0->'x', 1->'y', 2->'z'),
93// rotation axis index
94// and number of segments
95void halbach_array(const Vector &, Vector &);
96
97// A Field Boundary Condition for B = (Bx,By,Bz)
98static Vector b_uniform_(0);
99void a_bc_uniform(const Vector &, Vector&);
100
101// Phi_M Boundary Condition for H = (0,0,1)
103
104// Prints the program's logo to the given output stream
105void display_banner(ostream & os);
106
107int 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);
129 args.AddOption(&mesh_file, "-m", "--mesh",
130 "Mesh file to use.");
131 args.AddOption(&order, "-o", "--order",
132 "Finite element order (polynomial degree).");
133 args.AddOption(&serial_ref_levels, "-rs", "--serial-ref-levels",
134 "Number of serial refinement levels.");
135 args.AddOption(&parallel_ref_levels, "-rp", "--parallel-ref-levels",
136 "Number of parallel refinement levels.");
137 args.AddOption(&b_uniform_, "-ubbc", "--uniform-b-bc",
138 "Specify if the three components of the constant magnetic flux density");
139 args.AddOption(&pw_mu_, "-pwm", "--piecewise-mu",
140 "Piecewise values of Permeability");
141 args.AddOption(&ms_params_, "-ms", "--magnetic-shell-params",
142 "Center, Inner Radius, Outer Radius, and Permeability of Magnetic Shell");
143 args.AddOption(&cr_params_, "-cr", "--current-ring-params",
144 "Axis End Points, Inner Radius, Outer Radius and Total Current of Annulus");
145 args.AddOption(&bm_params_, "-bm", "--bar-magnet-params",
146 "Axis End Points, Radius, and Magnetic Field of Cylindrical Magnet");
147 args.AddOption(&ha_params_, "-ha", "--halbach-array-params",
148 "Bounding Box Corners and Number of Segments");
149 args.AddOption(&kbcs, "-kbcs", "--surface-current-bc",
150 "Surfaces for the Surface Current (K) Boundary Condition");
151 args.AddOption(&vbcs, "-vbcs", "--voltage-bc-surf",
152 "Voltage Boundary Condition Surfaces (to drive K)");
153 args.AddOption(&vbcv, "-vbcv", "--voltage-bc-vals",
154 "Voltage Boundary Condition Values (to drive K)");
155 args.AddOption(&maxit, "-maxit", "--max-amr-iterations",
156 "Max number of iterations in the main AMR loop.");
157 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
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 real_t local_max_err = errors.Max();
336 real_t global_max_err;
337 MPI_Allreduce(&local_max_err, &global_max_err, 1,
338 MPITypeMap<real_t>::mpi_type, MPI_MAX, pmesh.GetComm());
339
340 // Refine the elements whose error is larger than a fraction of the
341 // maximum element error.
342 const real_t frac = 0.5;
343 real_t 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
366void 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.
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_.
409{
410 real_t 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.
427void 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 real_t h = a.Norml2();
447
448 if ( h == 0.0 )
449 {
450 return;
451 }
452
453 real_t ra = cr_params_[2*x.Size()+0];
454 real_t rb = cr_params_[2*x.Size()+1];
455 if ( ra > rb )
456 {
457 real_t rc = ra;
458 ra = rb;
459 rb = rc;
460 }
461 real_t xa = xu*a;
462
463 if ( h > 0.0 )
464 {
465 xu.Add(-xa/(h*h),a);
466 }
467
468 real_t 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
477 j.Add(cr_params_[2*x.Size()+2]/(h*(rb-ra)),ju);
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.
484void 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 real_t h = a.Norml2();
501
502 if ( h == 0.0 )
503 {
504 return;
505 }
506
507 real_t r = bm_params_[2*x.Size()];
508 real_t xa = xu*a;
509
510 if ( h > 0.0 )
511 {
512 xu.Add(-xa/(h*h),a);
513 }
514
515 real_t xp = xu.Norml2();
516
517 if ( xa >= 0.0 && xa <= h*h && xp <= r )
518 {
519 m.Add(bm_params_[2*x.Size()+1]/h,a);
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.
526void 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 = static_cast<int>(n * (x[ai] - ha_params_[ai]) /
544 (ha_params_[ai+3] - ha_params_[ai]));
545
546 m[(ri + 1 + (i % 2)) % 3] = static_cast<int>(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).
551void 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).
562{
563 return -x(x.Size()-1);
564}
int Size() const
Return the logical size of the array.
Definition array.hpp:144
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
A coefficient that is constant across space and time.
A general function coefficient.
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Mesh data type.
Definition mesh.hpp:56
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition mesh.hpp:290
bool Nonconforming() const
Definition mesh.hpp:2229
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
bool RefineByError(const Array< real_t > &elem_error, real_t threshold, int nonconforming=-1, int nc_limit=0)
Definition mesh.cpp:10695
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition mesh.cpp:10626
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition mesh.cpp:6211
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
static int WorldSize()
Return the size of MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
A piecewise constant coefficient with the constants keyed off the element attribute numbers.
Class for parallel meshes.
Definition pmesh.hpp:34
MPI_Comm GetComm() const
Definition pmesh.hpp:402
void Rebalance()
Definition pmesh.cpp:4009
void Finalize(bool refine=false, bool fix_orientation=false) override
Finalize the construction of a general Mesh.
Definition pmesh.cpp:1524
Vector data type.
Definition vector.hpp:80
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:832
real_t Max() const
Returns the maximal element of the vector.
Definition vector.cpp:922
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
Vector & Add(const real_t a, const Vector &Va)
(*this) += a * Va
Definition vector.cpp:247
Data collection with VisIt I/O routines.
void RegisterVisItFields(VisItDataCollection &visit_dc)
int main()
real_t a
Definition lissajous.cpp:41
float real_t
Definition config.hpp:43
Helper struct to convert a C++ type to an MPI type.
real_t magnetic_shell(const Vector &)
Definition tesla.cpp:408
real_t magnetic_shell_inv(const Vector &x)
Definition tesla.cpp:78
void bar_magnet(const Vector &, Vector &)
Definition tesla.cpp:484
void display_banner(ostream &os)
Definition tesla.cpp:366
Coefficient * SetupInvPermeabilityCoefficient()
Definition tesla.cpp:379
void a_bc_uniform(const Vector &, Vector &)
Definition tesla.cpp:551
void halbach_array(const Vector &, Vector &)
Definition tesla.cpp:526
real_t phi_m_bc_uniform(const Vector &x)
Definition tesla.cpp:561
void current_ring(const Vector &, Vector &)
Definition tesla.cpp:427