MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
joule.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11//
12// -----------------------------------------------------
13// 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
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
112using namespace std;
113using namespace mfem;
114using namespace mfem::common;
115using namespace mfem::electromagnetics;
116
117void display_banner(ostream & os);
118
119static real_t mj_ = 0.0;
120static real_t sj_ = 0.0;
121static real_t wj_ = 0.0;
122
123// Initialize variables used in joule_solver.cpp
126
127int 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 real_t t_final = 100.0;
144 real_t dt = 0.5;
145 real_t mu = 1.0;
146 real_t sigma = 2.0*M_PI*10;
147 real_t Tcapacity = 1.0;
148 real_t Tconductivity = 0.01;
149 real_t 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);
161 "Mesh file to use.");
163 "Number of times to refine the mesh uniformly in serial.");
165 "Number of times to refine the mesh uniformly in parallel.");
167 "Order (degree) of the finite elements.");
169 "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3\n\t."
170 "\t 22 - Mid-Point, 23 - SDIRK23, 34 - SDIRK34.");
172 "Final time; start time is 0.");
174 "Time step.");
176 "Magnetic permeability coefficient.");
178 "Conductivity coefficient.");
180 "Frequency of oscillation.");
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.");
187 "Visualize every n-th timestep.");
189 "Name of the visit dump files");
191 "Print results (grid functions) to disk.");
193 "Enable AMR");
195 "Enable static condensation");
197 "Print matrices and vectors to disk");
199 "Hypre print level");
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, real_t> sigmaMap, InvTcondMap, TcapMap, InvTcapMap;
237 real_t sigmaAir = 0.0; // init to suppress gcc warning
238 real_t TcondAir;
239 real_t 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, real_t>(1, sigma));
256 sigmaMap.insert(pair<int, real_t>(2, sigmaAir));
257 sigmaMap.insert(pair<int, real_t>(3, sigmaAir));
258
259 InvTcondMap.insert(pair<int, real_t>(1, 1.0/Tconductivity));
260 InvTcondMap.insert(pair<int, real_t>(2, 1.0/TcondAir));
261 InvTcondMap.insert(pair<int, real_t>(3, 1.0/TcondAir));
262
263 TcapMap.insert(pair<int, real_t>(1, Tcapacity));
264 TcapMap.insert(pair<int, real_t>(2, TcapAir));
265 TcapMap.insert(pair<int, real_t>(3, TcapAir));
266
267 InvTcapMap.insert(pair<int, real_t>(1, 1.0/Tcapacity));
268 InvTcapMap.insert(pair<int, real_t>(2, 1.0/TcapAir));
269 InvTcapMap.insert(pair<int, real_t>(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.
442
443 ParFiniteElementSpace L2FESpace(pmesh, &L2FEC);
444 ParFiniteElementSpace HCurlFESpace(pmesh, &HCurlFEC);
445 ParFiniteElementSpace HDivFESpace(pmesh, &HDivFEC);
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();
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();
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]);
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
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,
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 real_t 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
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 real_t 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
727namespace mfem
728{
729
730namespace electromagnetics
731{
732
733void edot_bc(const Vector &x, Vector &E)
734{
735 E = 0.0;
736}
737
738void e_exact(const Vector &x, real_t t, Vector &E)
739{
740 E[0] = 0.0;
741 E[1] = 0.0;
742 E[2] = 0.0;
743}
744
745void b_exact(const Vector &x, real_t t, Vector &B)
746{
747 B[0] = 0.0;
748 B[1] = 0.0;
749 B[2] = 0.0;
750}
751
753{
754 real_t T = 0.0;
755 return T;
756}
757
759{
760 // the value
761 real_t 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
778void display_banner(ostream & os)
779{
780 os << " ____. .__ " << endl
781 << " | | ____ __ __| | ____ " << endl
782 << " | |/ _ \\| | \\ | _/ __ \\ " << endl
783 << "/\\__| ( <_> ) | / |_\\ ___/ " << endl
784 << "\\________|\\____/|____/|____/\\___ >" << endl
785 << " \\/ " << endl
786 << flush;
787}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
void DeleteAll()
Delete the whole array.
Definition array.hpp:864
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
Backward Euler ODE solver. L-stable.
Definition ode.hpp:412
A class to handle Vectors in a block fashion.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetTime(real_t t)
Set physical time (for time-dependent simulations)
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
A general function coefficient.
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Implicit midpoint method. A-stable, not L-stable.
Definition ode.hpp:425
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition mesh.cpp:10558
int GetAttribute(int i) const
Return the attribute of element i.
Definition mesh.hpp:1333
bool Nonconforming() const
Definition mesh.hpp:2229
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition mesh.cpp:10626
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition fe_coll.hpp:465
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition ode.hpp:24
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition ode.cpp:18
virtual void Step(Vector &x, real_t &t, real_t &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in].
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Abstract parallel finite element space.
Definition pfespace.hpp:29
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
Class for parallel grid function.
Definition pgridfunc.hpp:33
void Save(std::ostream &out) const override
void MakeRef(FiniteElementSpace *f, real_t *v) override
Make the ParGridFunction reference external data on a new FiniteElementSpace.
Class for parallel meshes.
Definition pmesh.hpp:34
MPI_Comm GetComm() const
Definition pmesh.hpp:402
void Rebalance()
Definition pmesh.cpp:4009
void Finalize(bool refine=false, bool fix_orientation=false) override
Finalize the construction of a general Mesh.
Definition pmesh.cpp:1524
void Print(std::ostream &out=mfem::out, const std::string &comments="") const override
Definition pmesh.cpp:4801
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:386
virtual void SetTime(real_t t)
Set the time for time dependent coefficients.
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
Data collection with VisIt I/O routines.
virtual void Save() override
Save the collection and a VisIt root file.
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.
void Debug(const char *basefilename, real_t time)
real_t ElectricLosses(ParGridFunction &E_gf) const
int close()
Close the socketstream.
real_t sigma(const Vector &x)
Definition maxwell.cpp:102
int problem
Definition ex15.cpp:62
int dim
Definition ex24.cpp:53
real_t freq
Definition ex24.cpp:54
real_t mu
Definition ex25.cpp:140
void E_exact(const Vector &, Vector &)
Definition ex3.cpp:251
int main()
HYPRE_Int HYPRE_BigInt
void display_banner(ostream &os)
Definition joule.cpp:778
int Ww
int Wy
int Wx
int offx
int Wh
int offy
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)
real_t p_bc(const Vector &x, real_t t)
Definition joule.cpp:758
void e_exact(const Vector &x, real_t t, Vector &E)
Definition joule.cpp:738
real_t t_exact(const Vector &x)
Definition joule.cpp:752
void edot_bc(const Vector &x, Vector &E)
Definition joule.cpp:733
void b_exact(const Vector &x, real_t t, Vector &B)
Definition joule.cpp:745
const int visport
void mfem_error(const char *msg)
Definition error.cpp:154
float real_t
Definition config.hpp:43
const char vishost[]
RefCoord t[3]