MFEM  v4.6.0
Finite element discretization library
joule.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 // Joule Miniapp: Transient Magnetics and Joule Heating
14 // -----------------------------------------------------
15 //
16 // This miniapp solves a time dependent eddy current problem, resulting in Joule
17 // heating.
18 //
19 // This version has electrostatic potential, Phi, which is a source term in the
20 // EM diffusion equation. The potential itself is driven by essential BC's
21 //
22 // Div sigma Grad Phi = 0
23 // sigma E = Curl B/mu - sigma grad Phi
24 // dB/dt = - Curl E
25 // F = -k Grad T
26 // c dT/dt = -Div(F) + sigma E.E,
27 //
28 // where B is the magnetic flux, E is the electric field, T is the temperature,
29 // F is the thermal flux, sigma is electrical conductivity, mu is the magnetic
30 // permeability, and alpha is the thermal diffusivity. The geometry of the
31 // domain is assumed to be as follows:
32 //
33 // boundary attribute 3
34 // +---------------------+
35 // boundary --->| | boundary
36 // attribute 1 | | attribute 2
37 // (driven) +---------------------+
38 //
39 // The voltage BC condition is essential BC on attribute 1 (front) and 2 (rear)
40 // given by function p_bc() at bottom of this file.
41 //
42 // The E-field boundary condition specifies the essential BC (n cross E) on
43 // attribute 1 (front) and 2 (rear) given by function edot_bc at bottom of this
44 // file. The E-field can be set on attribute 3 also.
45 //
46 // The thermal boundary condition for the flux F is the natural BC on attribute 1
47 // (front) and 2 (rear). This means that dT/dt = 0 on the boundaries, and the
48 // initial T = 0.
49 //
50 // See Section 3 for how the material propertied are assigned to mesh
51 // attributes, this needs to be changed for different applications.
52 //
53 // See Section 5 for how the boundary conditions are assigned to mesh
54 // attributes, this needs to be changed for different applications.
55 //
56 // This code supports a simple version of AMR, all elements containing material
57 // attribute 1 are (optionally) refined.
58 //
59 // Compile with: make joule
60 //
61 // Sample runs:
62 //
63 // mpirun -np 8 joule -m cylinder-hex.mesh -p rod
64 // mpirun -np 8 joule -m cylinder-tet.mesh -sc 1 -amr 1 -p rod
65 // mpirun -np 8 joule -m cylinder-hex-q2.gen -s 22 -dt 0.1 -tf 240.0 -p rod
66 //
67 // Options:
68 //
69 // -m [string] the mesh file name
70 // -o [int] the order of the basis
71 // -rs [int] number of times to serially refine the mesh
72 // -rp [int] number of times to refine the mesh in parallel
73 // -s [int] time integrator 1=Backward Euler, 2=SDIRK2, 3=SDIRK3,
74 // 22=Midpoint, 23=SDIRK23, 34=SDIRK34
75 // -tf [double] the final time
76 // -dt [double] time step
77 // -mu [double] the magnetic permeability
78 // -cnd [double] the electrical conductivity
79 // -f [double] the frequency of the applied EM BC
80 // -vis [int] GLVis -vis = true -no-vis = false
81 // -vs [int] visualization step
82 // -k [string] base file name for output file
83 // -print [int] print solution (gridfunctions) to disk 0 = no, 1 = yes
84 // -amr [int] 0 = no amr, 1 = amr
85 // -sc [int] 0 = no static condensation, 1 = use static condensation
86 // -p [string] specify the problem to run, "rod", "coil", or "test"
87 //
88 //
89 // NOTE: Example meshes for this miniapp are the included cylinder/rod meshes:
90 // cylinder-hex.mesh, cylinder-tet.mesh, cylinder-hex-q2.gen,
91 // cylinder-tet-q2.gen, as well as the coil.gen mesh which can be
92 // downloaded from github.com/mfem/data (its size is 21MB). Note that the
93 // meshes with the "gen" extension require MFEM to be built with NetCDF.
94 //
95 // NOTE: This code is set up to solve two example problems, 1) a straight metal
96 // rod surrounded by air, 2) a metal rod surrounded by a metal coil all
97 // surrounded by air. To specify problem (1) use the command line options
98 // "-p rod -m cylinder-hex-q2.gen", to specify problem (2) use the
99 // command line options "-p coil -m coil.gen". Problem (1) has two
100 // materials and problem (2) has three materials, and the BC's are
101 // different.
102 //
103 // NOTE: We write out, optionally, grid functions for P, E, B, W, F, and
104 // T. These can be visualized using "glvis -np 4 -m mesh.mesh -g E",
105 // assuming we used 4 processors.
106 
107 #include "joule_solver.hpp"
108 #include <memory>
109 #include <iostream>
110 #include <fstream>
111 
112 using namespace std;
113 using namespace mfem;
114 using namespace mfem::common;
115 using namespace mfem::electromagnetics;
116 
117 void display_banner(ostream & os);
118 
119 static double mj_ = 0.0;
120 static double sj_ = 0.0;
121 static double wj_ = 0.0;
122 
123 // Initialize variables used in joule_solver.cpp
126 
127 int main(int argc, char *argv[])
128 {
129  // 1. Initialize MPI and HYPRE.
130  Mpi::Init(argc, argv);
131  int myid = Mpi::WorldRank();
132  Hypre::Init();
133 
134  // print the cool banner
135  if (Mpi::Root()) { display_banner(cout); }
136 
137  // 2. Parse command-line options.
138  const char *mesh_file = "cylinder-hex.mesh";
139  int ser_ref_levels = 0;
140  int par_ref_levels = 0;
141  int order = 2;
142  int ode_solver_type = 1;
143  double t_final = 100.0;
144  double dt = 0.5;
145  double mu = 1.0;
146  double sigma = 2.0*M_PI*10;
147  double Tcapacity = 1.0;
148  double Tconductivity = 0.01;
149  double freq = 1.0/60.0;
150  bool visualization = true;
151  bool visit = true;
152  int vis_steps = 1;
153  int gfprint = 0;
154  const char *basename = "Joule";
155  int amr = 0;
156  int debug = 0;
157  const char *problem = "rod";
158 
159  OptionsParser args(argc, argv);
160  args.AddOption(&mesh_file, "-m", "--mesh",
161  "Mesh file to use.");
162  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
163  "Number of times to refine the mesh uniformly in serial.");
164  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
165  "Number of times to refine the mesh uniformly in parallel.");
166  args.AddOption(&order, "-o", "--order",
167  "Order (degree) of the finite elements.");
168  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
169  "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3\n\t."
170  "\t 22 - Mid-Point, 23 - SDIRK23, 34 - SDIRK34.");
171  args.AddOption(&t_final, "-tf", "--t-final",
172  "Final time; start time is 0.");
173  args.AddOption(&dt, "-dt", "--time-step",
174  "Time step.");
175  args.AddOption(&mu, "-mu", "--permeability",
176  "Magnetic permeability coefficient.");
177  args.AddOption(&sigma, "-cnd", "--sigma",
178  "Conductivity coefficient.");
179  args.AddOption(&freq, "-f", "--frequency",
180  "Frequency of oscillation.");
181  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
182  "--no-visualization",
183  "Enable or disable GLVis visualization.");
184  args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit",
185  "Enable or disable VisIt visualization.");
186  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
187  "Visualize every n-th timestep.");
188  args.AddOption(&basename, "-k", "--outputfilename",
189  "Name of the visit dump files");
190  args.AddOption(&gfprint, "-print", "--print",
191  "Print results (grid functions) to disk.");
192  args.AddOption(&amr, "-amr", "--amr",
193  "Enable AMR");
194  args.AddOption(&STATIC_COND, "-sc", "--static-condensation",
195  "Enable static condensation");
196  args.AddOption(&debug, "-debug", "--debug",
197  "Print matrices and vectors to disk");
198  args.AddOption(&SOLVER_PRINT_LEVEL, "-hl", "--hypre-print-level",
199  "Hypre print level");
200  args.AddOption(&problem, "-p", "--problem",
201  "Name of problem to run");
202  args.Parse();
203  if (!args.Good())
204  {
205  if (Mpi::Root())
206  {
207  args.PrintUsage(cout);
208  }
209  return 1;
210  }
211  if (Mpi::Root())
212  {
213  args.PrintOptions(cout);
214  }
215 
216  mj_ = mu;
217  sj_ = sigma;
218  wj_ = 2.0*M_PI*freq;
219 
220  if (Mpi::Root())
221  {
222  cout << "\nSkin depth sqrt(2.0/(wj*mj*sj)) = " << sqrt(2.0/(wj_*mj_*sj_))
223  << "\nSkin depth sqrt(2.0*dt/(mj*sj)) = " << sqrt(2.0*dt/(mj_*sj_))
224  << endl;
225  }
226 
227  // 3. Here material properties are assigned to mesh attributes. This code is
228  // not general, it is assumed the mesh has 3 regions each with a different
229  // integer attribute: 1, 2 or 3.
230  //
231  // The coil problem has three regions: 1) coil, 2) air, 3) the rod.
232  // The rod problem has two regions: 1) rod, 2) air.
233  //
234  // We can use the same material maps for both problems.
235 
236  std::map<int, double> sigmaMap, InvTcondMap, TcapMap, InvTcapMap;
237  double sigmaAir = 0.0; // init to suppress gcc warning
238  double TcondAir;
239  double TcapAir;
240  if (strcmp(problem,"rod")==0 || strcmp(problem,"coil")==0)
241  {
242  sigmaAir = 1.0e-6 * sigma;
243  TcondAir = 1.0e6 * Tconductivity;
244  TcapAir = 1.0 * Tcapacity;
245  }
246  else
247  {
248  cerr << "Problem " << problem << " not recognized\n";
249  mfem_error();
250  }
251 
252  if (strcmp(problem,"rod")==0 || strcmp(problem,"coil")==0)
253  {
254 
255  sigmaMap.insert(pair<int, double>(1, sigma));
256  sigmaMap.insert(pair<int, double>(2, sigmaAir));
257  sigmaMap.insert(pair<int, double>(3, sigmaAir));
258 
259  InvTcondMap.insert(pair<int, double>(1, 1.0/Tconductivity));
260  InvTcondMap.insert(pair<int, double>(2, 1.0/TcondAir));
261  InvTcondMap.insert(pair<int, double>(3, 1.0/TcondAir));
262 
263  TcapMap.insert(pair<int, double>(1, Tcapacity));
264  TcapMap.insert(pair<int, double>(2, TcapAir));
265  TcapMap.insert(pair<int, double>(3, TcapAir));
266 
267  InvTcapMap.insert(pair<int, double>(1, 1.0/Tcapacity));
268  InvTcapMap.insert(pair<int, double>(2, 1.0/TcapAir));
269  InvTcapMap.insert(pair<int, double>(3, 1.0/TcapAir));
270  }
271  else
272  {
273  cerr << "Problem " << problem << " not recognized\n";
274  mfem_error();
275  }
276 
277  // 4. Read the serial mesh from the given mesh file on all processors. We can
278  // handle triangular, quadrilateral, tetrahedral and hexahedral meshes
279  // with the same code.
280  Mesh *mesh;
281  mesh = new Mesh(mesh_file, 1, 1);
282  int dim = mesh->Dimension();
283 
284  // 5. Assign the boundary conditions
285  Array<int> ess_bdr(mesh->bdr_attributes.Max());
286  Array<int> thermal_ess_bdr(mesh->bdr_attributes.Max());
287  Array<int> poisson_ess_bdr(mesh->bdr_attributes.Max());
288  if (strcmp(problem,"coil")==0)
289  {
290  // BEGIN CODE FOR THE COIL PROBLEM
291  // For the coil in a box problem we have surfaces 1) coil end (+),
292  // 2) coil end (-), 3) five sides of box, 4) side of box with coil BC
293 
294  ess_bdr = 0;
295  ess_bdr[0] = 1; // boundary attribute 4 (index 3) is fixed
296  ess_bdr[1] = 1; // boundary attribute 4 (index 3) is fixed
297  ess_bdr[2] = 1; // boundary attribute 4 (index 3) is fixed
298  ess_bdr[3] = 1; // boundary attribute 4 (index 3) is fixed
299 
300  // Same as above, but this is for the thermal operator for HDiv
301  // formulation the essential BC is the flux
302 
303  thermal_ess_bdr = 0;
304  thermal_ess_bdr[2] = 1; // boundary attribute 4 (index 3) is fixed
305 
306  // Same as above, but this is for the poisson eq for H1 formulation the
307  // essential BC is the value of Phi
308 
309  poisson_ess_bdr = 0;
310  poisson_ess_bdr[0] = 1; // boundary attribute 1 (index 0) is fixed
311  poisson_ess_bdr[1] = 1; // boundary attribute 2 (index 1) is fixed
312  // END CODE FOR THE COIL PROBLEM
313  }
314  else if (strcmp(problem,"rod")==0)
315  {
316  // BEGIN CODE FOR THE STRAIGHT ROD PROBLEM
317  // the boundary conditions below are for the straight rod problem
318 
319  ess_bdr = 0;
320  ess_bdr[0] = 1; // boundary attribute 1 (index 0) is fixed (front)
321  ess_bdr[1] = 1; // boundary attribute 2 (index 1) is fixed (rear)
322  ess_bdr[2] = 1; // boundary attribute 3 (index 2) is fixed (outer)
323 
324  // Same as above, but this is for the thermal operator. For HDiv
325  // formulation the essential BC is the flux, which is zero on the front
326  // and sides. Note the Natural BC is T = 0 on the outer surface.
327 
328  thermal_ess_bdr = 0;
329  thermal_ess_bdr[0] = 1; // boundary attribute 1 (index 0) is fixed (front)
330  thermal_ess_bdr[1] = 1; // boundary attribute 2 (index 1) is fixed (rear)
331 
332  // Same as above, but this is for the poisson eq for H1 formulation the
333  // essential BC is the value of Phi
334 
335  poisson_ess_bdr = 0;
336  poisson_ess_bdr[0] = 1; // boundary attribute 1 (index 0) is fixed (front)
337  poisson_ess_bdr[1] = 1; // boundary attribute 2 (index 1) is fixed (back)
338  // END CODE FOR THE STRAIGHT ROD PROBLEM
339  }
340  else
341  {
342  cerr << "Problem " << problem << " not recognized\n";
343  mfem_error();
344  }
345 
346  // The following is required for mesh refinement
347  mesh->EnsureNCMesh();
348 
349  // 6. Define the ODE solver used for time integration. Several implicit
350  // methods are available, including singly diagonal implicit Runge-Kutta
351  // (SDIRK).
352  ODESolver *ode_solver;
353  switch (ode_solver_type)
354  {
355  // Implicit L-stable methods
356  case 1: ode_solver = new BackwardEulerSolver; break;
357  case 2: ode_solver = new SDIRK23Solver(2); break;
358  case 3: ode_solver = new SDIRK33Solver; break;
359  // Implicit A-stable methods (not L-stable)
360  case 22: ode_solver = new ImplicitMidpointSolver; break;
361  case 23: ode_solver = new SDIRK23Solver; break;
362  case 34: ode_solver = new SDIRK34Solver; break;
363  default:
364  if (Mpi::Root())
365  {
366  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
367  }
368  delete mesh;
369  return 3;
370  }
371 
372  // 7. Refine the mesh in serial to increase the resolution. In this example
373  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
374  // a command-line parameter.
375  for (int lev = 0; lev < ser_ref_levels; lev++)
376  {
377  mesh->UniformRefinement();
378  }
379 
380  // 8. Define a parallel mesh by a partitioning of the serial mesh. Refine
381  // this mesh further in parallel to increase the resolution. Once the
382  // parallel mesh is defined, the serial mesh can be deleted.
383  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
384  delete mesh;
385  for (int lev = 0; lev < par_ref_levels; lev++)
386  {
387  pmesh->UniformRefinement();
388  }
389  // Make sure tet-only meshes are marked for local refinement.
390  pmesh->Finalize(true);
391 
392  // 9. Apply non-uniform non-conforming mesh refinement to the mesh. The
393  // whole metal region is refined once, before the start of the time loop,
394  // i.e. this is not based on any error estimator.
395  if (amr == 1)
396  {
397  Array<int> ref_list;
398  int numElems = pmesh->GetNE();
399  for (int ielem = 0; ielem < numElems; ielem++)
400  {
401  int thisAtt = pmesh->GetAttribute(ielem);
402  if (thisAtt == 1)
403  {
404  ref_list.Append(ielem);
405  }
406  }
407 
408  pmesh->GeneralRefinement(ref_list);
409 
410  ref_list.DeleteAll();
411  }
412 
413  // 10. Rebalance the mesh. Since the mesh was adaptively refined in a
414  // non-uniform way it will be computationally unbalanced.
415  if (pmesh->Nonconforming())
416  {
417  pmesh->Rebalance();
418  }
419 
420  // 11. Define the parallel finite element spaces. We use:
421  //
422  // H(curl) for electric field,
423  // H(div) for magnetic flux,
424  // H(div) for thermal flux,
425  // H(grad)/H1 for electrostatic potential,
426  // L2 for temperature
427 
428  // L2 contains discontinuous "cell-center" finite elements, type 2 is
429  // "positive"
430  L2_FECollection L2FEC(order-1, dim);
431 
432  // ND contains Nedelec "edge-centered" vector finite elements with continuous
433  // tangential component.
434  ND_FECollection HCurlFEC(order, dim);
435 
436  // RT contains Raviart-Thomas "face-centered" vector finite elements with
437  // continuous normal component.
438  RT_FECollection HDivFEC(order-1, dim);
439 
440  // H1 contains continuous "node-centered" Lagrange finite elements.
441  H1_FECollection HGradFEC(order, dim);
442 
443  ParFiniteElementSpace L2FESpace(pmesh, &L2FEC);
444  ParFiniteElementSpace HCurlFESpace(pmesh, &HCurlFEC);
445  ParFiniteElementSpace HDivFESpace(pmesh, &HDivFEC);
446  ParFiniteElementSpace HGradFESpace(pmesh, &HGradFEC);
447 
448  // The terminology is TrueVSize is the unique (non-redundant) number of dofs
449  HYPRE_BigInt glob_size_l2 = L2FESpace.GlobalTrueVSize();
450  HYPRE_BigInt glob_size_nd = HCurlFESpace.GlobalTrueVSize();
451  HYPRE_BigInt glob_size_rt = HDivFESpace.GlobalTrueVSize();
452  HYPRE_BigInt glob_size_h1 = HGradFESpace.GlobalTrueVSize();
453 
454  if (Mpi::Root())
455  {
456  cout << "Number of Temperature Flux unknowns: " << glob_size_rt << endl;
457  cout << "Number of Temperature unknowns: " << glob_size_l2 << endl;
458  cout << "Number of Electric Field unknowns: " << glob_size_nd << endl;
459  cout << "Number of Magnetic Field unknowns: " << glob_size_rt << endl;
460  cout << "Number of Electrostatic unknowns: " << glob_size_h1 << endl;
461  }
462 
463  int Vsize_l2 = L2FESpace.GetVSize();
464  int Vsize_nd = HCurlFESpace.GetVSize();
465  int Vsize_rt = HDivFESpace.GetVSize();
466  int Vsize_h1 = HGradFESpace.GetVSize();
467 
468  // 12. Declare storage for field data.
469  // The big BlockVector stores the fields as
470  // 0 Temperature
471  // 1 Temperature Flux
472  // 2 P field
473  // 3 E field
474  // 4 B field
475  // 5 Joule Heating
476 
477  Array<int> true_offset(7);
478  true_offset[0] = 0;
479  true_offset[1] = true_offset[0] + Vsize_l2;
480  true_offset[2] = true_offset[1] + Vsize_rt;
481  true_offset[3] = true_offset[2] + Vsize_h1;
482  true_offset[4] = true_offset[3] + Vsize_nd;
483  true_offset[5] = true_offset[4] + Vsize_rt;
484  true_offset[6] = true_offset[5] + Vsize_l2;
485 
486  // The BlockVector is a large contiguous chunk of memory for storing required
487  // data for the hypre vectors, in this case: the temperature L2, the T-flux
488  // HDiv, the E-field HCurl, and the B-field HDiv, and scalar potential P.
489  BlockVector F(true_offset);
490 
491  // grid functions E, B, T, F, P, and w which is the Joule heating
492  ParGridFunction E_gf, B_gf, T_gf, F_gf, w_gf, P_gf;
493  T_gf.MakeRef(&L2FESpace,F, true_offset[0]);
494  F_gf.MakeRef(&HDivFESpace,F, true_offset[1]);
495  P_gf.MakeRef(&HGradFESpace,F,true_offset[2]);
496  E_gf.MakeRef(&HCurlFESpace,F,true_offset[3]);
497  B_gf.MakeRef(&HDivFESpace,F, true_offset[4]);
498  w_gf.MakeRef(&L2FESpace,F, true_offset[5]);
499 
500  // 13. Get the boundary conditions, set up the exact solution grid functions
501  // These VectorCoefficients have an Eval function. Note that e_exact and
502  // b_exact in this case are exact analytical solutions, taking a 3-vector
503  // point as input and returning a 3-vector field
506  FunctionCoefficient T_exact(t_exact);
507  E_exact.SetTime(0.0);
508  B_exact.SetTime(0.0);
509 
510  // 14. Initialize the Diffusion operator, the GLVis visualization and print
511  // the initial energies.
512  MagneticDiffusionEOperator oper(true_offset[6], L2FESpace, HCurlFESpace,
513  HDivFESpace, HGradFESpace,
514  ess_bdr, thermal_ess_bdr, poisson_ess_bdr,
515  mu, sigmaMap, TcapMap, InvTcapMap,
516  InvTcondMap);
517 
518  // This function initializes all the fields to zero or some provided IC
519  oper.Init(F);
520 
521  socketstream vis_T, vis_E, vis_B, vis_w, vis_P;
522  char vishost[] = "localhost";
523  int visport = 19916;
524  if (visualization)
525  {
526  // Make sure all ranks have sent their 'v' solution before initiating
527  // another set of GLVis connections (one from each rank):
528  MPI_Barrier(pmesh->GetComm());
529 
530  vis_T.precision(8);
531  vis_E.precision(8);
532  vis_B.precision(8);
533  vis_P.precision(8);
534  vis_w.precision(8);
535 
536  int Wx = 0, Wy = 0; // window position
537  int Ww = 350, Wh = 350; // window size
538  int offx = Ww+10, offy = Wh+45; // window offsets
539 
541  P_gf, "Electric Potential (Phi)", Wx, Wy, Ww, Wh);
542  Wx += offx;
543 
545  E_gf, "Electric Field (E)", Wx, Wy, Ww, Wh);
546  Wx += offx;
547 
549  B_gf, "Magnetic Field (B)", Wx, Wy, Ww, Wh);
550  Wx = 0;
551  Wy += offy;
552 
554  w_gf, "Joule Heating", Wx, Wy, Ww, Wh);
555 
556  Wx += offx;
557 
559  T_gf, "Temperature", Wx, Wy, Ww, Wh);
560  }
561  // VisIt visualization
562  VisItDataCollection visit_dc(basename, pmesh);
563  if ( visit )
564  {
565  visit_dc.RegisterField("E", &E_gf);
566  visit_dc.RegisterField("B", &B_gf);
567  visit_dc.RegisterField("T", &T_gf);
568  visit_dc.RegisterField("w", &w_gf);
569  visit_dc.RegisterField("Phi", &P_gf);
570  visit_dc.RegisterField("F", &F_gf);
571 
572  visit_dc.SetCycle(0);
573  visit_dc.SetTime(0.0);
574  visit_dc.Save();
575  }
576 
577  E_exact.SetTime(0.0);
578  B_exact.SetTime(0.0);
579 
580  // 15. Perform time-integration (looping over the time iterations, ti, with a
581  // time-step dt). The object oper is the MagneticDiffusionOperator which
582  // has a Mult() method and an ImplicitSolve() method which are used by
583  // the time integrators.
584  ode_solver->Init(oper);
585  double t = 0.0;
586 
587  bool last_step = false;
588  for (int ti = 1; !last_step; ti++)
589  {
590  if (t + dt >= t_final - dt/2)
591  {
592  last_step = true;
593  }
594 
595  // F is the vector of dofs, t is the current time, and dt is the time step
596  // to advance.
597  ode_solver->Step(F, t, dt);
598 
599  if (debug == 1)
600  {
601  oper.Debug(basename,t);
602  }
603 
604  if (gfprint == 1)
605  {
606  ostringstream T_name, E_name, B_name, F_name, w_name, P_name, mesh_name;
607  T_name << basename << "_" << setfill('0') << setw(6) << t << "_"
608  << "T." << setfill('0') << setw(6) << myid;
609  E_name << basename << "_" << setfill('0') << setw(6) << t << "_"
610  << "E." << setfill('0') << setw(6) << myid;
611  B_name << basename << "_" << setfill('0') << setw(6) << t << "_"
612  << "B." << setfill('0') << setw(6) << myid;
613  F_name << basename << "_" << setfill('0') << setw(6) << t << "_"
614  << "F." << setfill('0') << setw(6) << myid;
615  w_name << basename << "_" << setfill('0') << setw(6) << t << "_"
616  << "w." << setfill('0') << setw(6) << myid;
617  P_name << basename << "_" << setfill('0') << setw(6) << t << "_"
618  << "P." << setfill('0') << setw(6) << myid;
619  mesh_name << basename << "_" << setfill('0') << setw(6) << t << "_"
620  << "mesh." << setfill('0') << setw(6) << myid;
621 
622  ofstream mesh_ofs(mesh_name.str().c_str());
623  mesh_ofs.precision(8);
624  pmesh->Print(mesh_ofs);
625  mesh_ofs.close();
626 
627  ofstream T_ofs(T_name.str().c_str());
628  T_ofs.precision(8);
629  T_gf.Save(T_ofs);
630  T_ofs.close();
631 
632  ofstream E_ofs(E_name.str().c_str());
633  E_ofs.precision(8);
634  E_gf.Save(E_ofs);
635  E_ofs.close();
636 
637  ofstream B_ofs(B_name.str().c_str());
638  B_ofs.precision(8);
639  B_gf.Save(B_ofs);
640  B_ofs.close();
641 
642  ofstream F_ofs(F_name.str().c_str());
643  F_ofs.precision(8);
644  F_gf.Save(B_ofs);
645  F_ofs.close();
646 
647  ofstream P_ofs(P_name.str().c_str());
648  P_ofs.precision(8);
649  P_gf.Save(P_ofs);
650  P_ofs.close();
651 
652  ofstream w_ofs(w_name.str().c_str());
653  w_ofs.precision(8);
654  w_gf.Save(w_ofs);
655  w_ofs.close();
656  }
657 
658  if (last_step || (ti % vis_steps) == 0)
659  {
660  double el = oper.ElectricLosses(E_gf);
661 
662  if (Mpi::Root())
663  {
664  cout << fixed;
665  cout << "step " << setw(6) << ti << ",\tt = " << setw(6)
666  << setprecision(3) << t
667  << ",\tdot(E, J) = " << setprecision(8) << el << endl;
668  }
669 
670  // Make sure all ranks have sent their 'v' solution before initiating
671  // another set of GLVis connections (one from each rank):
672  MPI_Barrier(pmesh->GetComm());
673 
674  if (visualization)
675  {
676  int Wx = 0, Wy = 0; // window position
677  int Ww = 350, Wh = 350; // window size
678  int offx = Ww+10, offy = Wh+45; // window offsets
679 
681  P_gf, "Electric Potential (Phi)", Wx, Wy, Ww, Wh);
682  Wx += offx;
683 
685  E_gf, "Electric Field (E)", Wx, Wy, Ww, Wh);
686  Wx += offx;
687 
689  B_gf, "Magnetic Field (B)", Wx, Wy, Ww, Wh);
690 
691  Wx = 0;
692  Wy += offy;
693 
695  w_gf, "Joule Heating", Wx, Wy, Ww, Wh);
696 
697  Wx += offx;
698 
700  T_gf, "Temperature", Wx, Wy, Ww, Wh);
701  }
702 
703  if (visit)
704  {
705  visit_dc.SetCycle(ti);
706  visit_dc.SetTime(t);
707  visit_dc.Save();
708  }
709  }
710  }
711  if (visualization)
712  {
713  vis_T.close();
714  vis_E.close();
715  vis_B.close();
716  vis_w.close();
717  vis_P.close();
718  }
719 
720  // 16. Free the used memory.
721  delete ode_solver;
722  delete pmesh;
723 
724  return 0;
725 }
726 
727 namespace mfem
728 {
729 
730 namespace electromagnetics
731 {
732 
733 void edot_bc(const Vector &x, Vector &E)
734 {
735  E = 0.0;
736 }
737 
738 void e_exact(const Vector &x, double t, Vector &E)
739 {
740  E[0] = 0.0;
741  E[1] = 0.0;
742  E[2] = 0.0;
743 }
744 
745 void b_exact(const Vector &x, double t, Vector &B)
746 {
747  B[0] = 0.0;
748  B[1] = 0.0;
749  B[2] = 0.0;
750 }
751 
752 double t_exact(const Vector &x)
753 {
754  double T = 0.0;
755  return T;
756 }
757 
758 double p_bc(const Vector &x, double t)
759 {
760  // the value
761  double T;
762  if (x[2] < 0.0)
763  {
764  T = 1.0;
765  }
766  else
767  {
768  T = -1.0;
769  }
770 
771  return T*cos(wj_ * t);
772 }
773 
774 } // namespace electromagnetics
775 
776 } // namespace mfem
777 
778 void display_banner(ostream & os)
779 {
780  os << " ____. .__ " << endl
781  << " | | ____ __ __| | ____ " << endl
782  << " | |/ _ \\| | \\ | _/ __ \\ " << endl
783  << "/\\__| ( <_> ) | / |_\\ ___/ " << endl
784  << "\\________|\\____/|____/|____/\\___ >" << endl
785  << " \\/ " << endl
786  << flush;
787 }
int visport
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
A class to handle Vectors in a block fashion.
Definition: blockvector.hpp:30
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
int offx
virtual void Step(Vector &x, double &t, double &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
void e_exact(const Vector &x, double t, Vector &E)
Definition: joule.cpp:738
int Wx
bool Nonconforming() const
Definition: mesh.hpp:1969
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:23
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
Abstract parallel finite element space.
Definition: pfespace.hpp:28
STL namespace.
double ElectricLosses(ParGridFunction &E_gf) const
Backward Euler ODE solver. L-stable.
Definition: ode.hpp:411
void DeleteAll()
Delete the whole array.
Definition: array.hpp:854
void Debug(const char *basefilename, double time)
int close()
Close the socketstream.
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the ParGridFunction reference external data on a new FiniteElementSpace.
Definition: pgridfunc.cpp:111
int GetAttribute(int i) const
Return the attribute of element i.
Definition: mesh.hpp:1190
double p_bc(const Vector &x, double t)
Definition: joule.cpp:758
void Rebalance()
Definition: pmesh.cpp:4044
void Finalize(bool refine=false, bool fix_orientation=false) override
Finalize the construction of a general Mesh.
Definition: pmesh.cpp:1516
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
char vishost[]
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:759
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition: mesh.cpp:9888
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
Data collection with VisIt I/O routines.
void SetTime(double t)
Set physical time (for time-dependent simulations)
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:281
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
int main(int argc, char *argv[])
Definition: joule.cpp:127
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:380
A general vector function coefficient.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
virtual void Save() override
Save the collection and a VisIt root file.
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
int problem
Definition: ex15.cpp:62
HYPRE_Int HYPRE_BigInt
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:909
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
Implicit midpoint method. A-stable, not L-stable.
Definition: ode.hpp:424
virtual void SetTime(double t)
Set the time for time dependent coefficients.
int Wh
double mu
Definition: ex25.cpp:140
int dim
Definition: ex24.cpp:53
double freq
Definition: ex24.cpp:54
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:709
void E_exact(const Vector &, Vector &)
Definition: ex3.cpp:242
void edot_bc(const Vector &x, Vector &E)
Definition: joule.cpp:733
RefCoord t[3]
A general function coefficient.
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:454
Vector data type.
Definition: vector.hpp:58
double t_exact(const Vector &x)
Definition: joule.cpp:752
int offy
int Wy
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
double sigma(const Vector &x)
Definition: maxwell.cpp:102
int Ww
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4825
void b_exact(const Vector &x, double t, Vector &B)
Definition: joule.cpp:745
Class for parallel grid function.
Definition: pgridfunc.hpp:32
void display_banner(ostream &os)
Definition: joule.cpp:778
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:9820
Class for parallel meshes.
Definition: pmesh.hpp:32
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.