39#define jn(n, x) _jn(n, x)
40#define yn(n, x) _yn(n, x)
81 Array<int> * GetMarkedPMLElements() {
return &elems;}
84 void SetAttributes(
ParMesh *pmesh);
87 void StretchFunction(
const Vector &x, vector<complex<real_t>> &dxs);
97 PMLDiagMatrixCoefficient(
int dim,
void(*F)(
const Vector &, PML *,
112 (*Function)(transip, pml, K);
145template <
typename T> T
pow2(
const T &x) {
return x*x; }
157int main(
int argc,
char *argv[])
166 const char *mesh_file =
nullptr;
169 int par_ref_levels = 2;
172 bool herm_conv =
true;
173 bool slu_solver =
false;
174 bool mumps_solver =
false;
175 bool strumpack_solver =
false;
176 bool visualization = 1;
178 const char *device_config =
"cpu";
181 args.
AddOption(&mesh_file,
"-m",
"--mesh",
182 "Mesh file to use.");
184 "Finite element order (polynomial degree).");
185 args.
AddOption(&iprob,
"-prob",
"--problem",
"Problem case"
186 " 0: beam, 1: disc, 2: lshape, 3: fichera, 4: General");
187 args.
AddOption(&ref_levels,
"-rs",
"--refinements-serial",
188 "Number of serial refinements");
189 args.
AddOption(&par_ref_levels,
"-rp",
"--refinements-parallel",
190 "Number of parallel refinements");
192 "Permeability of free space (or 1/(spring constant)).");
194 "Permittivity of free space (or mass constant).");
196 "Frequency (in Hz).");
197 args.
AddOption(&herm_conv,
"-herm",
"--hermitian",
"-no-herm",
198 "--no-hermitian",
"Use convention for Hermitian operators.");
199#ifdef MFEM_USE_SUPERLU
200 args.
AddOption(&slu_solver,
"-slu",
"--superlu",
"-no-slu",
201 "--no-superlu",
"Use the SuperLU Solver.");
204 args.
AddOption(&mumps_solver,
"-mumps",
"--mumps-solver",
"-no-mumps",
205 "--no-mumps-solver",
"Use the MUMPS Solver.");
207#ifdef MFEM_USE_STRUMPACK
208 args.
AddOption(&strumpack_solver,
"-strumpack",
"--strumpack-solver",
209 "-no-strumpack",
"--no-strumpack-solver",
210 "Use the STRUMPACK Solver.");
212 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
213 "--no-visualization",
214 "Enable or disable GLVis visualization.");
215 args.
AddOption(&pa,
"-pa",
"--partial-assembly",
"-no-pa",
216 "--no-partial-assembly",
"Enable Partial Assembly.");
217 args.
AddOption(&device_config,
"-d",
"--device",
218 "Device configuration string, see Device::Configure().");
220 if (slu_solver + mumps_solver + strumpack_solver > 1)
223 cout <<
"WARNING: More than one of SuperLU, MUMPS, and STRUMPACK have"
224 <<
" been selected, please choose only one." << endl
225 <<
" Defaulting to SuperLU." << endl;
226 mumps_solver =
false;
227 strumpack_solver =
false;
230 if (iprob > 4) { iprob = 4; }
235 Device device(device_config);
236 if (myid == 0) { device.
Print(); }
245 mesh_file =
"../data/beam-hex.mesh";
248 mesh_file =
"../data/square-disc.mesh";
251 mesh_file =
"../data/l-shape.mesh";
254 mesh_file =
"../data/fichera.mesh";
258 mesh_file =
"../data/inline-quad.mesh";
276 Mesh * mesh =
new Mesh(mesh_file, 1, 1);
307 PML * pml =
new PML(mesh,length);
312 for (
int l = 0; l < ref_levels; l++)
321 for (
int l = 0; l < par_ref_levels; l++)
328 pml->SetAttributes(pmesh);
337 cout <<
"Number of finite element unknowns: " << size << endl;
352 for (
int j = 0; j < pmesh->
GetNBE(); j++)
362 if (center[0] == 1_r || center[0] == 0.5_r ||
369 if (center[0] == -1_r || center[0] == 0_r ||
370 center[1] == 0_r || center[2] == 0_r)
395 b.Vector::operator=(0.0);
421 attr = 0; attr[0] = 1;
439 int cdim = (
dim == 2) ? 1 :
dim;
464 if (pa) {
a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
469 a.FormLinearSystem(ess_tdof_list, x,
b, Ah, X, B);
472#ifdef MFEM_USE_SUPERLU
473 if (!pa && slu_solver)
487#ifdef MFEM_USE_STRUMPACK
488 if (!pa && strumpack_solver)
497 strumpack.
SetMatching(strumpack::MatchingJob::NONE);
501 strumpack.
Mult(B, X);
506 if (!pa && mumps_solver)
524 if (pa || (!slu_solver && !mumps_solver && !strumpack_solver))
555 std::unique_ptr<Operator> pc_r;
556 std::unique_ptr<Operator> pc_i;
591 a.RecoverFEMSolution(X,
b, x);
598 int order_quad = max(2, 2 * order + 1);
606 pml->GetMarkedPMLElements());
608 pml->GetMarkedPMLElements());
612 real_t norm_E_Re, norm_E_Im;
614 pml->GetMarkedPMLElements());
616 pml->GetMarkedPMLElements());
620 cout <<
"\n Relative Error (Re part): || E_h - E || / ||E|| = "
621 << L2Error_Re / norm_E_Re
622 <<
"\n Relative Error (Im part): || E_h - E || / ||E|| = "
623 << L2Error_Im / norm_E_Im
624 <<
"\n Total Error: "
625 << sqrt(L2Error_Re*L2Error_Re + L2Error_Im*L2Error_Im) <<
"\n\n";
632 ostringstream mesh_name, sol_r_name, sol_i_name;
633 mesh_name <<
"mesh." << setfill(
'0') << setw(6) << myid;
634 sol_r_name <<
"ex25p-sol_r." << setfill(
'0') << setw(6) << myid;
635 sol_i_name <<
"ex25p-sol_i." << setfill(
'0') << setw(6) << myid;
637 ofstream mesh_ofs(mesh_name.str().c_str());
638 mesh_ofs.precision(8);
639 pmesh->
Print(mesh_ofs);
641 ofstream sol_r_ofs(sol_r_name.str().c_str());
642 ofstream sol_i_ofs(sol_i_name.str().c_str());
643 sol_r_ofs.precision(8);
644 sol_i_ofs.precision(8);
654 keys = (
dim == 3) ?
"keys macF\n" : keys =
"keys amrRljcUUuu\n";
655 if (
prob ==
beam &&
dim == 3) {keys =
"keys macFFiYYYYYYYYYYYYYYYYYY\n";}
656 if (
prob ==
beam &&
dim == 2) {keys =
"keys amrRljcUUuuu\n"; }
663 sol_sock_re.precision(8);
664 sol_sock_re <<
"parallel " << num_procs <<
" " << myid <<
"\n"
665 <<
"solution\n" << *pmesh << x.
real() << keys
666 <<
"window_title 'Solution real part'" << flush;
667 MPI_Barrier(MPI_COMM_WORLD);
672 sol_sock_im.precision(8);
673 sol_sock_im <<
"parallel " << num_procs <<
" " << myid <<
"\n"
674 <<
"solution\n" << *pmesh << x.
imag() << keys
675 <<
"window_title 'Solution imag part'" << flush;
676 MPI_Barrier(MPI_COMM_WORLD);
684 sol_sock.precision(8);
685 sol_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
686 <<
"solution\n" << *pmesh << x_t << keys <<
"autoscale off\n"
687 <<
"window_title 'Harmonic Solution (t = 0.0 T)'"
688 <<
"pause\n" << flush;
692 cout <<
"GLVis visualization paused."
693 <<
" Press space (in the GLVis window) to resume it.\n";
702 oss <<
"Harmonic Solution (t = " <<
t <<
" T)";
706 sol_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n";
707 sol_sock <<
"solution\n" << *pmesh << x_t
708 <<
"window_title '" << oss.str() <<
"'" << flush;
726 for (
int i = 0; i <
dim; ++i)
729 r +=
pow2(x[i] - center[i]);
735 f[0] = coeff * exp(
alpha);
741 for (
int i = 0; i <
dim; ++i)
746 constexpr complex<real_t> zi = complex<real_t>(0., 1.);
762 real_t x0 = x(0) + shift(0);
763 real_t x1 = x(1) + shift(1);
764 real_t r = sqrt(x0 * x0 + x1 * x1);
768 complex<real_t> Ho, Ho_r, Ho_rr;
771 Ho_rr = -k * k * (1_r /
beta *
778 real_t r_xy = -(r_x / r) * r_y;
779 real_t r_xx = (1_r / r) * (1_r - r_x * r_x);
781 complex<real_t> val, val_xx, val_xy;
782 val =
real_t(0.25) * zi * Ho;
783 val_xx =
real_t(0.25) * zi * (r_xx * Ho_r + r_x * r_x * Ho_rr);
784 val_xy =
real_t(0.25) * zi * (r_xy * Ho_r + r_x * r_y * Ho_rr);
785 E[0] = zi / k * (k * k * val + val_xx);
786 E[1] = zi / k * val_xy;
790 real_t x0 = x(0) + shift(0);
791 real_t x1 = x(1) + shift(1);
792 real_t x2 = x(2) + shift(2);
793 real_t r = sqrt(x0 * x0 + x1 * x1 + x2 * x2);
798 real_t r_xx = (1_r / r) * (1_r - r_x * r_x);
799 real_t r_yx = -(r_y / r) * r_x;
800 real_t r_zx = -(r_z / r) * r_x;
802 complex<real_t> val, val_r, val_rr;
803 val = exp(zi * k * r) / r;
804 val_r = val / r * (zi * k * r - 1_r);
805 val_rr = val / (r * r) * (-k * k * r * r
808 complex<real_t> val_xx, val_yx, val_zx;
809 val_xx = val_rr * r_x * r_x + val_r * r_xx;
810 val_yx = val_rr * r_x * r_y + val_r * r_yx;
811 val_zx = val_rr * r_x * r_z + val_r * r_zx;
814 E[0] =
alpha * (k * k * val + val_xx);
815 E[1] =
alpha * val_yx;
816 E[2] =
alpha * val_zx;
826 E[1] = -zi * k / (
real_t) M_PI *
827 sin((
real_t) M_PI*x(2))*exp(zi * k10 * x(0));
831 E[1] = -zi * k / (
real_t) M_PI * exp(zi * k * x(0));
842 vector<complex<real_t>> Eval(E.
Size());
844 for (
int i = 0; i <
dim; ++i)
846 E[i] = Eval[i].real();
852 vector<complex<real_t>> Eval(E.
Size());
854 for (
int i = 0; i <
dim; ++i)
856 E[i] = Eval[i].imag();
865 for (
int i = 0; i <
dim; ++i)
877 vector<complex<real_t>> Eval(E.
Size());
879 for (
int i = 0; i <
dim; ++i)
881 E[i] = Eval[i].real();
892 for (
int i = 0; i <
dim; ++i)
904 vector<complex<real_t>> Eval(E.
Size());
906 for (
int i = 0; i <
dim; ++i)
908 E[i] = Eval[i].imag();
915 vector<complex<real_t>> dxs(
dim);
916 complex<real_t> det(1.0, 0.0);
917 pml->StretchFunction(x, dxs);
919 for (
int i = 0; i <
dim; ++i)
924 for (
int i = 0; i <
dim; ++i)
926 D(i) = (det /
pow2(dxs[i])).real();
932 vector<complex<real_t>> dxs(
dim);
933 complex<real_t> det = 1.0;
934 pml->StretchFunction(x, dxs);
936 for (
int i = 0; i <
dim; ++i)
941 for (
int i = 0; i <
dim; ++i)
943 D(i) = (det /
pow2(dxs[i])).imag();
949 vector<complex<real_t>> dxs(
dim);
950 complex<real_t> det = 1.0;
951 pml->StretchFunction(x, dxs);
953 for (
int i = 0; i <
dim; ++i)
958 for (
int i = 0; i <
dim; ++i)
960 D(i) = abs(det /
pow2(dxs[i]));
966 vector<complex<real_t>> dxs(
dim);
967 complex<real_t> det(1.0, 0.0);
968 pml->StretchFunction(x, dxs);
970 for (
int i = 0; i <
dim; ++i)
978 D = (1_r / det).real();
982 for (
int i = 0; i <
dim; ++i)
984 D(i) = (
pow2(dxs[i]) / det).real();
991 vector<complex<real_t>> dxs(
dim);
992 complex<real_t> det = 1.0;
993 pml->StretchFunction(x, dxs);
995 for (
int i = 0; i <
dim; ++i)
1002 D = (1_r / det).imag();
1006 for (
int i = 0; i <
dim; ++i)
1008 D(i) = (
pow2(dxs[i]) / det).imag();
1015 vector<complex<real_t>> dxs(
dim);
1016 complex<real_t> det = 1.0;
1017 pml->StretchFunction(x, dxs);
1019 for (
int i = 0; i <
dim; ++i)
1030 for (
int i = 0; i <
dim; ++i)
1032 D(i) = abs(
pow2(dxs[i]) / det);
1038 : mesh(mesh_), length(length_)
1040 dim = mesh->Dimension();
1044void PML::SetBoundaries()
1050 for (
int i = 0; i <
dim; i++)
1052 dom_bdr(i, 0) = pmin(i);
1053 dom_bdr(i, 1) = pmax(i);
1054 comp_dom_bdr(i, 0) = dom_bdr(i, 0) + length(i, 0);
1055 comp_dom_bdr(i, 1) = dom_bdr(i, 1) - length(i, 1);
1059void PML::SetAttributes(
ParMesh *pmesh)
1062 for (
int i = 0; i < pmesh->
GetNBE(); ++i)
1067 int nrelem = pmesh->
GetNE();
1073 for (
int i = 0; i < nrelem; ++i)
1076 bool in_pml =
false;
1083 int nrvert = vertices.
Size();
1086 for (
int iv = 0; iv < nrvert; ++iv)
1088 int vert_idx = vertices[iv];
1090 for (
int comp = 0; comp <
dim; ++comp)
1092 if (coords[comp] > comp_dom_bdr(comp, 1) ||
1093 coords[comp] < comp_dom_bdr(comp, 0))
1109void PML::StretchFunction(
const Vector &x,
1110 vector<complex<real_t>> &dxs)
1112 constexpr complex<real_t> zi = complex<real_t>(0., 1.);
1120 for (
int i = 0; i <
dim; ++i)
1125 coeff = n * c / k / pow(length(i, 1), n);
1126 dxs[i] = 1_r + zi * coeff *
1131 coeff = n * c / k / pow(length(i, 0), n);
1132 dxs[i] = 1_r + zi * coeff *
Dynamic 2D array using row-major layout.
void SetSize(int m, int n)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
int Size() const
Return the logical size of the array.
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
A class to handle Block diagonal preconditioners in a matrix-free implementation.
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
Specialization of the ComplexOperator built from a pair of HypreParMatrices.
@ HERMITIAN
Native convention for Hermitian operators.
@ BLOCK_SYMMETRIC
Alternate convention for damping operators.
A coefficient that is constant across space and time.
Integrator for for Nedelec elements.
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Abstract data type element.
virtual void GetVertices(Array< int > &v) const =0
Get the indices defining the vertices.
void SetAttribute(const int attr)
Set element's attribute.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the GMRES method.
const IntegrationPoint & GetCenter(int GeomType) const
Return the center of the given Geometry::Type, GeomType.
The Auxiliary-space Maxwell Solver in hypre.
Wrapper for hypre's ParCSR matrix class.
MPI_Comm GetComm() const
MPI communicator.
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Class for integration point with weight.
Class for an integration rule - an Array of IntegrationPoint.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
void SetRelTol(real_t rtol)
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
void SetMaxIter(int max_it)
void SetAbsTol(real_t atol)
MUMPS: A Parallel Sparse Direct Solver.
void SetPrintLevel(int print_lvl)
Set the error print level for MUMPS.
void Mult(const Vector &x, Vector &y) const
Solve .
void SetMatrixSymType(MatType mtype)
Set the matrix type.
void SetOperator(const Operator &op)
Set the Operator and perform factorization.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Geometry::Type GetBdrElementGeometry(int i) const
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
const Element * GetElement(int i) const
Return pointer to the i'th element object.
int GetNE() const
Returns number of elements.
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
int Dimension() const
Dimension of the reference space used within the elements.
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
ElementTransformation * GetBdrElementTransformation(int i)
Returns a pointer to the transformation defining the i-th boundary element.
int GetNBE() const
Returns number of boundary elements.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Array< int > attributes
A list of all unique element attributes used by the Mesh.
const real_t * GetVertex(int i) const
Return pointer to vertex i's coordinates.
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static int WorldSize()
Return the size of 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.
Pointer to an Operator of a specified type.
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get().
Jacobi smoothing for a given bilinear form (no matrix necessary).
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.
virtual void ProjectBdrCoefficientTangent(VectorCoefficient &real_coeff, VectorCoefficient &imag_coeff, Array< int > &attr)
Abstract parallel finite element space.
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
HYPRE_BigInt GlobalTrueVSize() const
int GetTrueVSize() const override
Return the number of local vector true dofs.
Class for parallel grid function.
void Save(std::ostream &out) const override
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
Class for parallel meshes.
void SetAttributes() override
Determine the sets of unique attribute values in domain and boundary elements.
void Print(std::ostream &out=mfem::out, const std::string &comments="") const override
Derived coefficient that takes the value of the parent coefficient for the active attributes and is z...
void SetOperator(const Operator &op)
Set the operator/matrix.
void SetMatching(strumpack::MatchingJob job)
Configure static pivoting for stability.
void SetCompression(strumpack::CompressionType type)
Select compression for sparse data types.
void SetPrintFactorStatistics(bool print_stat)
Set up verbose printing during the factor step.
void SetKrylovSolver(strumpack::KrylovSolver method)
Set the Krylov solver method to use.
void Mult(const Vector &x, Vector &y) const
Factor and solve the linear system .
void SetFromCommandLine()
Set options that were captured from the command line.
void SetPrintSolveStatistics(bool print_stat)
Set up verbose printing during the solve step.
void SetReorderingStrategy(strumpack::ReorderingStrategy method)
Set matrix reordering strategy.
Vector coefficient defined as a product of scalar and vector coefficients.
Scaled Operator B: x -> a A(x).
void SetColumnPermutation(superlu::ColPerm col_perm)
Specify how to permute the columns of the matrix.
void Mult(const Vector &x, Vector &y) const
Factor and solve the linear system .
void SetSymmetricPattern(bool sym)
Specify whether the matrix has a symmetric pattern to avoid extra work (default false)
void SetOperator(const Operator &op)
Set the operator/matrix.
void SetPrintStatistics(bool print_stat)
Specify whether to print the solver statistics (default true)
Base class for vector Coefficients that optionally depend on time and space.
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
for VectorFiniteElements (Nedelec, Raviart-Thomas)
A general vector function coefficient.
Derived vector coefficient that has the value of the parent vector where it is active and is zero oth...
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
void detJ_JT_J_inv_abs(const Vector &x, PML *pml, Vector &D)
void detJ_inv_JT_J_Im(const Vector &x, PML *pml, Vector &D)
void E_bdr_data_Re(const Vector &x, Vector &E)
void detJ_JT_J_inv_Re(const Vector &x, PML *pml, Vector &D)
void detJ_inv_JT_J_abs(const Vector &x, PML *pml, Vector &D)
void maxwell_solution(const Vector &x, vector< complex< real_t > > &Eval)
void source(const Vector &x, Vector &f)
Array2D< real_t > comp_domain_bdr
void E_exact_Re(const Vector &x, Vector &E)
void detJ_JT_J_inv_Im(const Vector &x, PML *pml, Vector &D)
void detJ_inv_JT_J_Re(const Vector &x, PML *pml, Vector &D)
void E_exact_Im(const Vector &x, Vector &E)
Array2D< real_t > domain_bdr
void E_bdr_data_Im(const Vector &x, Vector &E)
@ PARMETIS
Sequential ordering on structure of using the PARMETIS package.
void add(const Vector &v1, const Vector &v2, Vector &v)
std::function< real_t(const Vector &)> f(real_t mass_coeff)
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)