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