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