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