113 using namespace mfem;
114 using namespace mfem::electromagnetics;
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;
127 int main(
int argc,
char *argv[])
137 const char *mesh_file =
"cylinder-hex.mesh";
138 int ser_ref_levels = 0;
139 int par_ref_levels = 0;
141 int ode_solver_type = 1;
142 double t_final = 100.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;
155 const char *basename =
"Joule";
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.");
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",
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.");
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");
221 kj_ = sqrt(0.5*wj_*mj_*sj_);
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_)));
242 std::map<int, double> sigmaMap, InvTcondMap, TcapMap, InvTcapMap;
246 if (strcmp(problem,
"rod")==0 || strcmp(problem,
"coil")==0)
248 sigmaAir = 1.0e-6 * sigma;
249 TcondAir = 1.0e6 * Tconductivity;
250 TcapAir = 1.0 * Tcapacity;
254 cerr <<
"Problem " << problem <<
" not recognized\n";
258 if (strcmp(problem,
"rod")==0 || strcmp(problem,
"coil")==0)
261 sigmaMap.insert(pair<int, double>(1, sigma));
262 sigmaMap.insert(pair<int, double>(2, sigmaAir));
263 sigmaMap.insert(pair<int, double>(3, sigmaAir));
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));
269 TcapMap.insert(pair<int, double>(1, Tcapacity));
270 TcapMap.insert(pair<int, double>(2, TcapAir));
271 TcapMap.insert(pair<int, double>(3, TcapAir));
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));
279 cerr <<
"Problem " << problem <<
" not recognized\n";
287 mesh =
new Mesh(mesh_file, 1, 1);
294 if (strcmp(problem,
"coil")==0)
310 thermal_ess_bdr[2] = 1;
316 poisson_ess_bdr[0] = 1;
317 poisson_ess_bdr[1] = 1;
320 else if (strcmp(problem,
"rod")==0)
335 thermal_ess_bdr[0] = 1;
336 thermal_ess_bdr[1] = 1;
342 poisson_ess_bdr[0] = 1;
343 poisson_ess_bdr[1] = 1;
348 cerr <<
"Problem " << problem <<
" not recognized\n";
359 switch (ode_solver_type)
372 cout <<
"Unknown ODE solver type: " << ode_solver_type <<
'\n';
381 for (
int lev = 0; lev < ser_ref_levels; lev++)
391 for (
int lev = 0; lev < par_ref_levels; lev++)
402 int numElems = pmesh->
GetNE();
403 for (
int ielem = 0; ielem < numElems; ielem++)
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;
471 int Vsize_l2 = L2FESpace.
GetVSize();
472 int Vsize_nd = HCurlFESpace.
GetVSize();
473 int Vsize_rt = HDivFESpace.
GetVSize();
474 int Vsize_h1 = HGradFESpace.
GetVSize();
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;
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]);
520 HDivFESpace, HGradFESpace,
521 ess_bdr, thermal_ess_bdr, poisson_ess_bdr,
522 mu, sigmaMap, TcapMap, InvTcapMap,
529 char vishost[] =
"localhost";
544 int Ww = 350, Wh = 350;
545 int offx = Ww+10, offy = Wh+45;
548 P_gf,
"Electric Potential (Phi)", Wx, Wy, Ww, Wh);
552 E_gf,
"Electric Field (E)", Wx, Wy, Ww, Wh);
556 B_gf,
"Magnetic Field (B)", Wx, Wy, Ww, Wh);
561 w_gf,
"Joule Heating", Wx, Wy, Ww, Wh);
566 T_gf,
"Temperature", Wx, Wy, Ww, Wh);
591 ode_solver->
Init(oper);
594 bool last_step =
false;
595 for (
int ti = 1; !last_step; ti++)
597 if (t + dt >= t_final - dt/2)
604 ode_solver->
Step(F, t, dt);
608 oper.
Debug(basename,t);
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;
629 ofstream mesh_ofs(mesh_name.str().c_str());
630 mesh_ofs.precision(8);
631 pmesh->
Print(mesh_ofs);
634 ofstream T_ofs(T_name.str().c_str());
639 ofstream E_ofs(E_name.str().c_str());
644 ofstream B_ofs(B_name.str().c_str());
649 ofstream F_ofs(F_name.str().c_str());
654 ofstream P_ofs(P_name.str().c_str());
659 ofstream w_ofs(w_name.str().c_str());
665 if (last_step || (ti % vis_steps) == 0)
672 cout <<
"step " << setw(6) << ti <<
",\tt = " << setw(6)
673 << setprecision(3) << t
674 <<
",\tdot(E, J) = " << setprecision(8) << el << endl;
684 int Ww = 350, Wh = 350;
685 int offx = Ww+10, offy = Wh+45;
688 P_gf,
"Electric Potential (Phi)", Wx, Wy, Ww, Wh);
692 E_gf,
"Electric Field (E)", Wx, Wy, Ww, Wh);
696 B_gf,
"Magnetic Field (B)", Wx, Wy, Ww, Wh);
702 w_gf,
"Joule Heating", Wx, Wy, Ww, Wh);
707 T_gf,
"Temperature", Wx, Wy, Ww, Wh);
737 namespace electromagnetics
778 return T*cos(wj_ * t);
787 os <<
" ____. .__ " << endl
788 <<
" | | ____ __ __| | ____ " << endl
789 <<
" | |/ _ \\| | \\ | _/ __ \\ " << endl
790 <<
"/\\__| ( <_> ) | / |_\\ ___/ " << endl
791 <<
"\\________|\\____/|____/|____/\\___ >" << endl
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
virtual void Init(TimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
virtual void ReorientTetMesh()
See the remarks for the serial version in mesh.hpp.
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)
double t_exact(Vector &x)
HYPRE_Int GlobalTrueVSize()
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
int GetNE() const
Returns number of elements.
virtual void Save()
Save the collection and a VisIt root file.
virtual void Save(std::ostream &out) const
Abstract parallel finite element space.
int main(int argc, char *argv[])
Backward Euler ODE solver. L-stable.
void DeleteAll()
Delete whole array.
void Debug(const char *basefilename, double time)
double p_bc(const Vector &x, double t)
void MakeRef(ParFiniteElementSpace *f, double *v)
Make the ParGridFunction reference external data on a new ParFiniteElementSpace.
void Rebalance()
Load balance the mesh. NC meshes only.
bool Nonconforming() const
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.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
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)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
bool Root() const
Return true if WorldRank() == 0.
void PrintUsage(std::ostream &out) const
void SetTime(double t)
Set physical time (for time-dependent simulations)
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
double ElectricLosses(ParGridFunction &E_gf) const
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
int WorldRank() const
Return MPI_COMM_WORLD's rank.
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)
void mfem_error(const char *msg)
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.
void PrintOptions(std::ostream &out) const
void E_exact(const Vector &, Vector &)
void EnsureNCMesh(bool triangles_nonconforming=false)
void edot_bc(const Vector &x, Vector &E)
class for C-function coefficient
Arbitrary order H(curl)-conforming Nedelec finite elements.
Arbitrary order H1-conforming (continuous) finite elements.
void b_exact(const Vector &x, double t, Vector &B)
Class for parallel grid function.
void display_banner(ostream &os)
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Class for parallel meshes.
int GetAttribute(int i) const
Return the attribute of element i.
Arbitrary order "L2-conforming" discontinuous finite elements.
virtual void Print(std::ostream &out=std::cout) const