87static Vector ds_params_(0);
93static Vector ms_params_(0);
99static Vector cs_params_(0);
105static Vector dp_params_(0);
120static real_t tScale_ = 1e-9;
127int main(
int argc,
char *argv[])
129#ifdef MFEM_USE_SINGLE
130 cout <<
"This miniapp is not supported in single precision.\n\n";
131 return MFEM_SKIP_RETURN_VALUE;
140 const char *mesh_file =
"../../data/ball-nurbs.mesh";
143 int serial_ref_levels = 0;
144 int parallel_ref_levels = 0;
146 bool visualization =
true;
158 args.
AddOption(&mesh_file,
"-m",
"--mesh",
159 "Mesh file to use.");
160 args.
AddOption(&sOrder,
"-so",
"--spatial-order",
161 "Finite element order (polynomial degree).");
162 args.
AddOption(&tOrder,
"-to",
"--temporal-order",
163 "Time integration order.");
164 args.
AddOption(&serial_ref_levels,
"-rs",
"--serial-ref-levels",
165 "Number of serial refinement levels.");
166 args.
AddOption(¶llel_ref_levels,
"-rp",
"--parallel-ref-levels",
167 "Number of parallel refinement levels.");
168 args.
AddOption(&dtsf,
"-sf",
"--dt-safety-factor",
169 "Used to reduce the time step below the upper bound.");
170 args.
AddOption(&ti,
"-ti",
"--initial-time",
171 "Beginning of time interval to simulate (ns).");
172 args.
AddOption(&tf,
"-tf",
"--final-time",
173 "End of time interval to simulate (ns).");
174 args.
AddOption(&ts,
"-ts",
"--snapshot-time",
175 "Time between snapshots (ns).");
176 args.
AddOption(&ds_params_,
"-ds",
"--dielectric-sphere-params",
177 "Center, Radius, and Permittivity of Dielectric Sphere");
178 args.
AddOption(&ms_params_,
"-ms",
"--magnetic-shell-params",
179 "Center, Inner Radius, Outer Radius, and Permeability "
180 "of Magnetic Shell");
181 args.
AddOption(&cs_params_,
"-cs",
"--conductive-sphere-params",
182 "Center, Radius, and Conductivity of Conductive Sphere");
183 args.
AddOption(&dp_params_,
"-dp",
"--dipole-pulse-params",
184 "Axis End Points, Radius, Amplitude, "
185 "Pulse Center (ns), Pulse Width (ns)");
186 args.
AddOption(&abcs,
"-abcs",
"--absorbing-bc-surf",
187 "Absorbing Boundary Condition Surfaces");
188 args.
AddOption(&dbcs,
"-dbcs",
"--dirichlet-bc-surf",
189 "Dirichlet Boundary Condition Surfaces");
190 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
191 "--no-visualization",
192 "Enable or disable GLVis visualization.");
193 args.
AddOption(&visit,
"-visit",
"--visit",
"-no-visit",
194 "--no-visualization",
195 "Enable or disable VisIt visualization.");
196 args.
AddOption(&visport,
"-p",
"--send-port",
"Socket for GLVis.");
215 ifstream imesh(mesh_file);
220 cerr <<
"\nCan not open mesh file: " << mesh_file <<
'\n' << endl;
224 mesh =
new Mesh(imesh, 1, 1);
231 if (serial_ref_levels > 0) { serial_ref_levels--; }
238 for (
int l = 0; l < serial_ref_levels; l++)
246 ParMesh pmesh(MPI_COMM_WORLD, *mesh);
250 for (
int l = 0; l < parallel_ref_levels; l++)
258 ( ms_params_.
Size() > 0 ) ?
muInv : NULL,
259 ( cs_params_.
Size() > 0 ) ?
sigma : NULL,
260 ( dp_params_.
Size() > 0 ) ?
j_src : NULL,
279 cout <<
"Energy(" << ti <<
"ns): " << energy <<
"J" << endl;
292 cout <<
"Maximum Time Step: " << dtmax / tScale_ <<
"ns" << endl;
299 cout <<
"Number of Time Steps: " << nsteps << endl;
300 cout <<
"Time Step Size: " << dt / tScale_ <<
"ns" << endl;
343 max(t + dt, ti + ts * it));
349 cout <<
"Energy(" << t/tScale_ <<
"ns): " << energy <<
"J" << endl;
376 os <<
" ___ ____ " << endl
377 <<
" / | / / __ __ " << endl
378 <<
" / |_/ _ /__ ___ _____ _ __ ____ | | | | " << endl
379 <<
" / \\__ \\ \\ \\/ /\\ \\/ \\/ // __ \\| | | | "
381 <<
" / /|_/ // __ \\_> < \\ /\\ ___/| |_| |__ " << endl
382 <<
"/___/ /_ /(____ /__/\\_ \\ \\/\\_/ \\___ >____/____/ " << endl
383 <<
" \\/ \\/ \\/ \\/ " << endl
393 for (
int i=0; i<x.
Size(); i++)
395 r2 += (x(i)-ds_params_(i))*(x(i)-ds_params_(i));
398 if ( sqrt(r2) <= ds_params_(x.
Size()) )
400 return ds_params_(x.
Size()+1) * epsilon0_;
412 for (
int i=0; i<x.
Size(); i++)
414 r2 += (x(i)-ms_params_(i))*(x(i)-ms_params_(i));
417 if ( sqrt(r2) >= ms_params_(x.
Size()) &&
418 sqrt(r2) <= ms_params_(x.
Size()+1) )
420 return mu0_*ms_params_(x.
Size()+2);
431 for (
int i=0; i<x.
Size(); i++)
433 r2 += (x(i)-cs_params_(i))*(x(i)-cs_params_(i));
436 if ( sqrt(r2) <= cs_params_(x.
Size()) )
438 return cs_params_(x.
Size()+1);
448 MFEM_ASSERT(x.
Size() == 3,
"current source requires 3D space.");
458 for (
int i=0; i<x.
Size(); i++)
460 xu[i] -= dp_params_[i];
461 v[i] = dp_params_[x.
Size()+i] - dp_params_[i];
475 real_t c = dp_params_[2*x.
Size()+3] * tScale_;
484 if ( xv >= 0.0 && xv <= h && xp <= r )
489 j *=
a * (t -
b) * exp(-0.5 * pow((t-
b)/c, 2.0)) / (c * c);
516 real_t dsteps = tmax/dtmax;
518 int nsteps =
static_cast<int>(pow(10,(
int)ceil(log10(dsteps))));
520 for (
int i=1; i<=5; i++)
522 int a = (int)ceil(log10(dsteps/pow(5.0,i)));
523 int nstepsi = (int)pow(5,i)*max(1,(
int)pow(10,
a));
525 nsteps = min(nsteps,nstepsi);
int Size() const
Return the logical size of the array.
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
NURBSExtension * NURBSext
Optional NURBS mesh extension.
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
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...
bool Good() const
Return true if the command line options were parsed successfully.
Class for parallel meshes.
virtual void Init(Operator &P, TimeDependentOperator &F)
virtual void Run(Vector &q, Vector &p, real_t &t, real_t &dt, real_t tf)
Variable order Symplectic Integration Algorithm (orders 1-4)
virtual void SetTime(const real_t t_)
Set the current time.
A general vector function coefficient.
real_t Norml2() const
Returns the l2 norm of the vector.
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
Vector & Add(const real_t a, const Vector &Va)
(*this) += a * Va
Data collection with VisIt I/O routines.
real_t GetMaximumTimeStep() const
void WriteVisItFields(int it=0)
void SetInitialEField(VectorCoefficient &EFieldCoef)
void RegisterVisItFields(VisItDataCollection &visit_dc)
void SetInitialBField(VectorCoefficient &BFieldCoef)
void DisplayToGLVis(int visport=19916)
real_t sigma(const Vector &x)
void dipole_pulse(const Vector &x, real_t t, Vector &j)
real_t magnetic_shell(const Vector &)
void EFieldFunc(const Vector &, Vector &)
real_t muInv(const Vector &x)
int SnapTimeStep(real_t tmax, real_t dtmax, real_t &dt)
void j_src(const Vector &x, real_t t, Vector &j)
void dEdtBCFunc(const Vector &x, real_t t, Vector &E)
real_t dielectric_sphere(const Vector &)
void display_banner(ostream &os)
void BFieldFunc(const Vector &, Vector &)
real_t conductive_sphere(const Vector &)