Compare commits
5
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
 Will Pazner
|
552d6857cb | ||
|
Will Pazner
|
a37d46e917 | ||
|
Will Pazner
|
4acdb072b6 | ||
|
Will Pazner
|
9dbb184537 | ||
|
Will Pazner
|
d67762a1c9 |
@@ -55,11 +55,8 @@ doc/warnings.log
|
||||
|
||||
examples/ex[0-9]
|
||||
examples/ex[0-9]p
|
||||
examples/ex[0-9]-orth
|
||||
examples/ex[0-9]p-orth
|
||||
examples/ex1[04-9]
|
||||
examples/ex1[0-9]p
|
||||
examples/ex1[0-9]p-cyl
|
||||
examples/ex2[0-9]
|
||||
examples/ex2[0-9]p
|
||||
examples/ex3[0-9]
|
||||
|
||||
@@ -52,7 +52,6 @@ if (MFEM_USE_MPI)
|
||||
list(APPEND ALL_EXE_SRCS
|
||||
ex0p.cpp
|
||||
ex1p.cpp
|
||||
ex1p-orth.cpp
|
||||
ex2p.cpp
|
||||
ex3p.cpp
|
||||
ex4p.cpp
|
||||
@@ -63,11 +62,8 @@ if (MFEM_USE_MPI)
|
||||
ex9p.cpp
|
||||
ex10p.cpp
|
||||
ex11p.cpp
|
||||
ex11p-cyl.cpp
|
||||
ex12p.cpp
|
||||
ex13p.cpp
|
||||
ex13p-cyl.cpp
|
||||
ex13p-cyl-3d.cpp
|
||||
ex14p.cpp
|
||||
ex15p.cpp
|
||||
ex16p.cpp
|
||||
|
||||
@@ -1,399 +0,0 @@
|
||||
// MFEM Example 11-cyl - Parallel Version
|
||||
//
|
||||
// Compile with: make ex11p-cyl
|
||||
//
|
||||
// Sample runs: mpirun -np 4 ex11p-cyl
|
||||
// mpirun -np 4 ex11p-cyl -o 2
|
||||
// mpirun -np 4 ex11p-cyl -o 2 -e 0
|
||||
//
|
||||
// Description: This example code demonstrates the use of MFEM to solve PDEs
|
||||
// on an axisymmetric domain. The eigenvalue problem:
|
||||
// -Delta u = lambda u
|
||||
// with homogeneous Dirichlet boundary conditions is solved on
|
||||
// a cylindrical domain by meshing only a rectangle in the
|
||||
// rho, z plane. In cylindrical coordinates the weak form of
|
||||
// the eigenvalue problem is given by:
|
||||
// (rho Grad(u), Grad(v)) = lambda (rho u, v)
|
||||
//
|
||||
// We compute the five lowest eigenmodes by discretizing
|
||||
// the Laplacian and Mass operators using a FE space of the
|
||||
// specified order and compare to the known values. Because the
|
||||
// eigenvalue spectrum of a domain is unique this provides a
|
||||
// reliable test that the axisymmetric domain is being faithfully
|
||||
// characterized.
|
||||
//
|
||||
// The example highlights the use of specialized coefficients
|
||||
// with existing operators to mimic axisymmetric domains. The
|
||||
// gradient of each eigenmode is also computed and displayed to
|
||||
// ilustrate that no special steps need to be taken to compute
|
||||
// gradients in this coordinate system.
|
||||
//
|
||||
// We recommend viewing Example 11 before viewing this example.
|
||||
|
||||
#include "mfem.hpp"
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
|
||||
using namespace std;
|
||||
using namespace mfem;
|
||||
|
||||
// Zeros of Bessel function J_0
|
||||
static double J0z[] = {2.40482555769577,
|
||||
5.52007811028631,
|
||||
8.65372791291101,
|
||||
11.7915344390143
|
||||
};
|
||||
|
||||
// Modes numbers in the rho and z directions for the first five eigenmodes
|
||||
static int mode_nums[] = {0, 1,
|
||||
1, 1,
|
||||
0, 2,
|
||||
1, 2,
|
||||
2, 1
|
||||
};
|
||||
|
||||
double rhoFunc(const Vector &x)
|
||||
{
|
||||
return x[0];
|
||||
}
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
// 1. Initialize MPI.
|
||||
int num_procs, myid;
|
||||
MPI_Init(&argc, &argv);
|
||||
MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
|
||||
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
|
||||
|
||||
// 2. Parse command-line options.
|
||||
int nr = 1;
|
||||
int nz = 1;
|
||||
int el_type_flag = 1;
|
||||
Element::Type el_type;
|
||||
int ser_ref_levels = 2;
|
||||
int par_ref_levels = 1;
|
||||
int order = 1;
|
||||
int nev = 5;
|
||||
int seed = 75;
|
||||
bool slu_solver = false;
|
||||
bool sp_solver = false;
|
||||
bool visualization = 1;
|
||||
|
||||
OptionsParser args(argc, argv);
|
||||
args.AddOption(&nz, "-nz", "--num-elements-z",
|
||||
"Number of elements in z-direction.");
|
||||
args.AddOption(&nr, "-nr", "--num-elements-rho",
|
||||
"Number of elements in radial direction.");
|
||||
args.AddOption(&el_type_flag, "-e", "--element-type",
|
||||
"Element type: 0 - Triangle, 1 - Quadrilateral.");
|
||||
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",
|
||||
"Finite element order (polynomial degree) or -1 for"
|
||||
" isoparametric space.");
|
||||
args.AddOption(&seed, "-s", "--seed",
|
||||
"Random seed used to initialize LOBPCG.");
|
||||
#ifdef MFEM_USE_SUPERLU
|
||||
args.AddOption(&slu_solver, "-slu", "--superlu", "-no-slu",
|
||||
"--no-superlu", "Use the SuperLU Solver.");
|
||||
#endif
|
||||
#ifdef MFEM_USE_STRUMPACK
|
||||
args.AddOption(&sp_solver, "-sp", "--strumpack", "-no-sp",
|
||||
"--no-strumpack", "Use the STRUMPACK Solver.");
|
||||
#endif
|
||||
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
|
||||
"--no-visualization",
|
||||
"Enable or disable GLVis visualization.");
|
||||
args.Parse();
|
||||
if (slu_solver && sp_solver)
|
||||
{
|
||||
if (myid == 0)
|
||||
cout << "WARNING: Both SuperLU and STRUMPACK have been selected,"
|
||||
<< " please choose either one." << endl
|
||||
<< " Defaulting to SuperLU." << endl;
|
||||
sp_solver = false;
|
||||
}
|
||||
// The command line options are also passed to the STRUMPACK
|
||||
// solver. So do not exit if some options are not recognized.
|
||||
if (!sp_solver)
|
||||
{
|
||||
if (!args.Good())
|
||||
{
|
||||
if (myid == 0)
|
||||
{
|
||||
args.PrintUsage(cout);
|
||||
}
|
||||
MPI_Finalize();
|
||||
return 1;
|
||||
}
|
||||
}
|
||||
if (myid == 0)
|
||||
{
|
||||
args.PrintOptions(cout);
|
||||
}
|
||||
|
||||
// The output mesh could be quadrilaterals or triangles
|
||||
el_type = (el_type_flag == 0) ? Element::TRIANGLE : Element::QUADRILATERAL;
|
||||
if (el_type != Element::TRIANGLE && el_type != Element::QUADRILATERAL)
|
||||
{
|
||||
cout << "Unsupported element type" << endl;
|
||||
exit(1);
|
||||
}
|
||||
|
||||
// 3. Prepare a rectangular mesh with the desired dimensions and element
|
||||
// type. Other 2D meshes could be used but then we couldn't check the
|
||||
// eigenvalues.
|
||||
Mesh *mesh = new Mesh(nr, nz, el_type);
|
||||
int dim = mesh->Dimension();
|
||||
|
||||
// 4. Refine the serial mesh on all processors to increase the resolution. In
|
||||
// this example we do 'ref_levels' of uniform refinement (2 by default, or
|
||||
// specified on the command line with -rs).
|
||||
for (int lev = 0; lev < ser_ref_levels; lev++)
|
||||
{
|
||||
mesh->UniformRefinement();
|
||||
}
|
||||
|
||||
// 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
|
||||
// this mesh further in parallel to increase the resolution (1 time by
|
||||
// default, or specified on the command line with -rp). 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();
|
||||
}
|
||||
|
||||
// 6. Define a parallel finite element space on the parallel mesh. Here we
|
||||
// use continuous Lagrange finite elements (H1) of the specified order.
|
||||
// We also create a Nedelec space to represent the gradients of the modes.
|
||||
H1_FECollection fec_h1(order, dim);
|
||||
ND_FECollection fec_nd(order, dim);
|
||||
ParFiniteElementSpace fespace_h1(pmesh, &fec_h1);
|
||||
ParFiniteElementSpace fespace_nd(pmesh, &fec_nd);
|
||||
HYPRE_Int size = fespace_h1.GlobalTrueVSize();
|
||||
if (myid == 0)
|
||||
{
|
||||
cout << "Number of unknowns: " << size << endl;
|
||||
}
|
||||
|
||||
// 7. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
|
||||
// element space. The first corresponds to the Laplacian operator -Delta,
|
||||
// while the second is a simple mass matrix needed on the right hand side
|
||||
// of the generalized eigenvalue problem below. The boundary conditions
|
||||
// are implemented by elimination with special values on the diagonal to
|
||||
// shift the Dirichlet eigenvalues out of the computational range. After
|
||||
// serial and parallel assembly we extract the corresponding parallel
|
||||
// matrices A and M.
|
||||
FunctionCoefficient rhoCoef(rhoFunc);
|
||||
Array<int> ess_bdr(pmesh->bdr_attributes.Max());
|
||||
ess_bdr = 1; // Homogeneous Dirichlet BCs everywhere except for
|
||||
ess_bdr[3] = 0; // attribute 4 which is on the axis of symmetry.
|
||||
|
||||
ParBilinearForm *a = new ParBilinearForm(&fespace_h1);
|
||||
a->AddDomainIntegrator(new DiffusionIntegrator(rhoCoef));
|
||||
a->Assemble();
|
||||
a->EliminateEssentialBCDiag(ess_bdr, 1.0);
|
||||
a->Finalize();
|
||||
|
||||
ParBilinearForm *m = new ParBilinearForm(&fespace_h1);
|
||||
m->AddDomainIntegrator(new MassIntegrator(rhoCoef));
|
||||
m->Assemble();
|
||||
// shift the eigenvalue corresponding to eliminated dofs to a large value
|
||||
m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
|
||||
m->Finalize();
|
||||
|
||||
HypreParMatrix *A = a->ParallelAssemble();
|
||||
HypreParMatrix *M = m->ParallelAssemble();
|
||||
|
||||
#if defined(MFEM_USE_SUPERLU) || defined(MFEM_USE_STRUMPACK)
|
||||
Operator * Arow = NULL;
|
||||
#ifdef MFEM_USE_SUPERLU
|
||||
if (slu_solver)
|
||||
{
|
||||
Arow = new SuperLURowLocMatrix(*A);
|
||||
}
|
||||
#endif
|
||||
#ifdef MFEM_USE_STRUMPACK
|
||||
if (sp_solver)
|
||||
{
|
||||
Arow = new STRUMPACKRowLocMatrix(*A);
|
||||
}
|
||||
#endif
|
||||
#endif
|
||||
|
||||
delete a;
|
||||
delete m;
|
||||
|
||||
// 8. Define and configure the LOBPCG eigensolver and the BoomerAMG
|
||||
// preconditioner for A to be used within the solver. Set the matrices
|
||||
// which define the generalized eigenproblem A x = lambda M x.
|
||||
Solver * precond = NULL;
|
||||
if (!slu_solver && !sp_solver)
|
||||
{
|
||||
HypreBoomerAMG * amg = new HypreBoomerAMG(*A);
|
||||
amg->SetPrintLevel(0);
|
||||
precond = amg;
|
||||
}
|
||||
else
|
||||
{
|
||||
#ifdef MFEM_USE_SUPERLU
|
||||
if (slu_solver)
|
||||
{
|
||||
SuperLUSolver * superlu = new SuperLUSolver(MPI_COMM_WORLD);
|
||||
superlu->SetPrintStatistics(false);
|
||||
superlu->SetSymmetricPattern(true);
|
||||
superlu->SetColumnPermutation(superlu::PARMETIS);
|
||||
superlu->SetOperator(*Arow);
|
||||
precond = superlu;
|
||||
}
|
||||
#endif
|
||||
#ifdef MFEM_USE_STRUMPACK
|
||||
if (sp_solver)
|
||||
{
|
||||
STRUMPACKSolver * strumpack = new STRUMPACKSolver(argc, argv, MPI_COMM_WORLD);
|
||||
strumpack->SetPrintFactorStatistics(true);
|
||||
strumpack->SetPrintSolveStatistics(false);
|
||||
strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
|
||||
strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
|
||||
strumpack->DisableMatching();
|
||||
strumpack->SetOperator(*Arow);
|
||||
strumpack->SetFromCommandLine();
|
||||
precond = strumpack;
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
|
||||
lobpcg->SetNumModes(nev);
|
||||
lobpcg->SetRandomSeed(seed);
|
||||
lobpcg->SetPreconditioner(*precond);
|
||||
lobpcg->SetMaxIter(200);
|
||||
lobpcg->SetTol(1e-8);
|
||||
lobpcg->SetPrecondUsageMode(1);
|
||||
lobpcg->SetPrintLevel(1);
|
||||
lobpcg->SetMassMatrix(*M);
|
||||
lobpcg->SetOperator(*A);
|
||||
|
||||
// 9. Compute the eigenmodes and extract the array of eigenvalues. Define a
|
||||
// parallel grid function to represent each of the eigenmodes returned by
|
||||
// the solver. Also define a discrete gradient operator.
|
||||
Array<double> eigenvalues;
|
||||
lobpcg->Solve();
|
||||
lobpcg->GetEigenvalues(eigenvalues);
|
||||
ParGridFunction x(&fespace_h1);
|
||||
ParGridFunction dx(&fespace_nd);
|
||||
|
||||
ParDiscreteLinearOperator grad(&fespace_h1, &fespace_nd);
|
||||
grad.AddDomainInterpolator(new GradientInterpolator());
|
||||
grad.Assemble();
|
||||
|
||||
if ( myid == 0 )
|
||||
{
|
||||
// Display the eigenvalues and their relative errors
|
||||
cout << "\nRelative error in eigenvalues:\n";
|
||||
for (int i=0; i<nev; i++)
|
||||
{
|
||||
double lambda =
|
||||
pow(J0z[mode_nums[2*i]], 2) +
|
||||
pow(M_PI * mode_nums[2*i+1], 2);
|
||||
cout << "Lambda " << i+1 << '/' << nev << " = " << eigenvalues[i]
|
||||
<< ", rel err = " << fabs(eigenvalues[i] - lambda) / lambda
|
||||
<< endl;
|
||||
}
|
||||
cout << endl;
|
||||
}
|
||||
|
||||
// 10. Save the refined mesh and the modes in parallel. This output can be
|
||||
// viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
|
||||
{
|
||||
ostringstream mesh_name, mode_name;
|
||||
mesh_name << "mesh." << setfill('0') << setw(6) << myid;
|
||||
|
||||
ofstream mesh_ofs(mesh_name.str().c_str());
|
||||
mesh_ofs.precision(8);
|
||||
pmesh->Print(mesh_ofs);
|
||||
|
||||
for (int i=0; i<nev; i++)
|
||||
{
|
||||
// convert eigenvector from HypreParVector to ParGridFunction
|
||||
x = lobpcg->GetEigenvector(i);
|
||||
|
||||
mode_name << "mode_" << setfill('0') << setw(2) << i << "."
|
||||
<< setfill('0') << setw(6) << myid;
|
||||
|
||||
ofstream mode_ofs(mode_name.str().c_str());
|
||||
mode_ofs.precision(8);
|
||||
x.Save(mode_ofs);
|
||||
mode_name.str("");
|
||||
}
|
||||
}
|
||||
|
||||
// 11. Send the solution by socket to a GLVis server.
|
||||
if (visualization)
|
||||
{
|
||||
char vishost[] = "localhost";
|
||||
int visport = 19916;
|
||||
socketstream mode_sock(vishost, visport);
|
||||
mode_sock.precision(8);
|
||||
socketstream grad_sock(vishost, visport);
|
||||
grad_sock.precision(8);
|
||||
|
||||
for (int i=0; i<nev; i++)
|
||||
{
|
||||
if ( myid == 0 )
|
||||
{
|
||||
cout << "Eigenmode " << i+1 << '/' << nev
|
||||
<< ", Lambda = " << eigenvalues[i] << endl;
|
||||
}
|
||||
|
||||
// convert eigenvector from HypreParVector to ParGridFunction
|
||||
x = lobpcg->GetEigenvector(i);
|
||||
|
||||
grad.Mult(x, dx);
|
||||
|
||||
mode_sock << "parallel " << num_procs << " " << myid << "\n"
|
||||
<< "solution\n" << *pmesh << x << flush
|
||||
<< "window_title 'Eigenmode " << i+1 << '/' << nev
|
||||
<< ", Lambda = " << eigenvalues[i] << "'" << endl;
|
||||
|
||||
grad_sock << "parallel " << num_procs << " " << myid << "\n"
|
||||
<< "solution\n" << *pmesh << dx << flush
|
||||
<< "window_title 'Grad of Eigenmode " << i+1 << '/' << nev
|
||||
<< ", Lambda = " << eigenvalues[i] << "'"
|
||||
<< "window_geometry 400 0 400 350" << endl;
|
||||
|
||||
char c;
|
||||
if (myid == 0)
|
||||
{
|
||||
cout << "press (q)uit or (c)ontinue --> " << flush;
|
||||
cin >> c;
|
||||
}
|
||||
MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
|
||||
|
||||
if (c != 'c')
|
||||
{
|
||||
break;
|
||||
}
|
||||
}
|
||||
mode_sock.close();
|
||||
}
|
||||
|
||||
// 12. Free the used memory.
|
||||
delete lobpcg;
|
||||
delete precond;
|
||||
delete M;
|
||||
delete A;
|
||||
#if defined(MFEM_USE_SUPERLU) || defined(MFEM_USE_STRUMPACK)
|
||||
delete Arow;
|
||||
#endif
|
||||
delete pmesh;
|
||||
|
||||
MPI_Finalize();
|
||||
|
||||
return 0;
|
||||
}
|
||||
@@ -1,563 +0,0 @@
|
||||
// MFEM Example 13-cyl - Parallel Version
|
||||
//
|
||||
// Compile with: make ex13p-cyl
|
||||
//
|
||||
// Sample runs: mpirun -np 4 ex13p-cyl
|
||||
// mpirun -np 4 ex13p-cyl -o 2
|
||||
// mpirun -np 4 ex13p-cyl -o 2 -e 0
|
||||
//
|
||||
// Description: This example code demonstrates the use of MFEM to solve the
|
||||
// Maxwell (electromagnetic) eigenvalue problem on an
|
||||
// axisymmetric domain. The eigenvalue problem:
|
||||
// curl curl E = lambda E
|
||||
// with homogeneous Dirichlet boundary conditions E x n = 0 is
|
||||
// solved on a cylindrical domain by meshing only a rectangle in
|
||||
// the rho, z plane. In cylindrical coordinates the weak form of
|
||||
// the eigenvalue problem is given by:
|
||||
// (rho Curl(u), Curl(v)) = lambda (rho u, v)
|
||||
//
|
||||
// We compute the eight lowest nonzero eigenmodes by discretizing
|
||||
// the curl curl operator using a Nedelec FE space of the
|
||||
// specified order and compare to the known values. Because the
|
||||
// eigenvalue spectrum of a domain is unique this provides a
|
||||
// reliable test that the axisymmetric domain is being faithfully
|
||||
// characterized.
|
||||
//
|
||||
// In two dimensions the curl curl operator, with isotropic
|
||||
// material coefficients, splits into two separate PDEs. The rho
|
||||
// and z components form a 2D vector field in the rho-z plane
|
||||
// which is discretized with the 2D Nedelec vector basis
|
||||
// functions. The Maxwell eigenvalue problem for the rho-z field
|
||||
// can be written with cartesian operators as:
|
||||
// curl (rho curl E_rz) = lambda rho E_rz
|
||||
//
|
||||
// The angular component can be discretized with the 2D H1 scalar
|
||||
// basis functions. The angular portion of the eigenvalue problem
|
||||
// can be written with cartesian operators as:
|
||||
// -div (rho grad E_phi) + (1/rho) E_phi = lambda rho E_phi
|
||||
//
|
||||
// The example highlights the use of specialized coefficients
|
||||
// with existing operators to mimic axisymmetric domains. The
|
||||
// curl of each eigenmode is also computed and displayed to
|
||||
// ilustrate that no special steps need to be taken to compute
|
||||
// the curl in this coordinate system.
|
||||
//
|
||||
// We recommend viewing examples 13 and 11-cyl before viewing this
|
||||
// example.
|
||||
|
||||
#include "mfem.hpp"
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
|
||||
using namespace std;
|
||||
using namespace mfem;
|
||||
|
||||
// Zeros of Bessel function J_0
|
||||
static double J0z[] = {2.404825557695773,
|
||||
5.520078110286312,
|
||||
8.653727912911013,
|
||||
11.79153443901428
|
||||
};
|
||||
|
||||
// Zeros of Bessel function J_1
|
||||
static double J1z[] = {3.831705970207512,
|
||||
7.015586669815622,
|
||||
10.17346813506272,
|
||||
13.32369193631422
|
||||
};
|
||||
|
||||
// Modes numbers in the rho and z directions for the first five eigenmodes
|
||||
static int mode_nums_rz[] = {0, 0,
|
||||
0, 1,
|
||||
1, 0,
|
||||
1, 1,
|
||||
0, 2,
|
||||
1, 2,
|
||||
2, 0,
|
||||
2, 1
|
||||
};
|
||||
|
||||
// Modes numbers in the phi direction for the first five eigenmodes
|
||||
static int mode_nums_phi[] = {0, 1,
|
||||
0, 2,
|
||||
1, 1,
|
||||
1, 2,
|
||||
0, 3,
|
||||
2, 1,
|
||||
1, 3,
|
||||
2, 2
|
||||
};
|
||||
|
||||
// Mode polarizations for the first several modes; 0 - rz, 1 - phi
|
||||
static int mode_type[] = {0,0,1,0,0,0,1,1,0,0,0,1,0,1,1,0,1,0,1};
|
||||
|
||||
double rhoFunc(const Vector &x)
|
||||
{
|
||||
return x[0];
|
||||
}
|
||||
|
||||
double rhoInvFunc(const Vector &x)
|
||||
{
|
||||
return 1.0 / x[0];
|
||||
}
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
// 1. Initialize MPI.
|
||||
int num_procs, myid;
|
||||
MPI_Init(&argc, &argv);
|
||||
MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
|
||||
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
|
||||
|
||||
// 2. Parse command-line options.
|
||||
int dim = 2;
|
||||
int nr = 1;
|
||||
int nz = 1;
|
||||
int el_type_flag = 1;
|
||||
Element::Type el_type;
|
||||
int ser_ref_levels = 2;
|
||||
int par_ref_levels = 1;
|
||||
int order = 1;
|
||||
int nev = 8;
|
||||
bool visualization = 1;
|
||||
|
||||
OptionsParser args(argc, argv);
|
||||
args.AddOption(&nz, "-nz", "--num-elements-z",
|
||||
"Number of elements in z-direction.");
|
||||
args.AddOption(&nr, "-nr", "--num-elements-rho",
|
||||
"Number of elements in radial direction.");
|
||||
args.AddOption(&el_type_flag, "-e", "--element-type",
|
||||
"Element type: 0 - Triangle, 1 - Quadrilateral.");
|
||||
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",
|
||||
"Finite element order (polynomial degree) or -1 for"
|
||||
" isoparametric space.");
|
||||
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
|
||||
"--no-visualization",
|
||||
"Enable or disable GLVis visualization.");
|
||||
args.Parse();
|
||||
if (!args.Good())
|
||||
{
|
||||
if (myid == 0)
|
||||
{
|
||||
args.PrintUsage(cout);
|
||||
}
|
||||
MPI_Finalize();
|
||||
return 1;
|
||||
}
|
||||
if (myid == 0)
|
||||
{
|
||||
args.PrintOptions(cout);
|
||||
}
|
||||
|
||||
// The output mesh could be quadrilaterals or triangles
|
||||
el_type = (el_type_flag == 0) ? Element::TRIANGLE : Element::QUADRILATERAL;
|
||||
if (el_type != Element::TRIANGLE && el_type != Element::QUADRILATERAL)
|
||||
{
|
||||
cout << "Unsupported element type" << endl;
|
||||
exit(1);
|
||||
}
|
||||
|
||||
// 3. Prepare a rectangular mesh with the desired dimensions and element
|
||||
// type. Other 2D meshes could be used but then we couldn't check the
|
||||
// eigenvalues.
|
||||
ParMesh pmesh;
|
||||
{
|
||||
Mesh mesh = Mesh::MakeCartesian2D(nr, nz, el_type);
|
||||
|
||||
// 4. Refine the serial mesh on all processors to increase the resolution. In
|
||||
// this example we do 'ref_levels' of uniform refinement (2 by default, or
|
||||
// specified on the command line with -rs).
|
||||
for (int lev = 0; lev < ser_ref_levels; lev++)
|
||||
{
|
||||
mesh.UniformRefinement();
|
||||
}
|
||||
|
||||
// 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
|
||||
// this mesh further in parallel to increase the resolution (1 time by
|
||||
// default, or specified on the command line with -rp). Once the parallel
|
||||
// mesh is defined, the serial mesh can be deleted.
|
||||
pmesh = ParMesh(MPI_COMM_WORLD, mesh);
|
||||
for (int lev = 0; lev < par_ref_levels; lev++)
|
||||
{
|
||||
pmesh.UniformRefinement();
|
||||
}
|
||||
pmesh.ReorientTetMesh();
|
||||
}
|
||||
|
||||
// 6. Define a parallel finite element space on the parallel mesh. Here we
|
||||
// use the Nedelec finite elements (ND) of the specified order. We also
|
||||
// create an L2 space to represent the z-component of the curl of the
|
||||
// modes.
|
||||
ND_FECollection fec_nd(order, dim);
|
||||
RT_FECollection fec_rt(order - 1, dim);
|
||||
H1_FECollection fec_ndp(order, dim);
|
||||
L2_FECollection fec_rtp(order - 1, dim,
|
||||
BasisType::GaussLegendre, FiniteElement::INTEGRAL);
|
||||
L2_FECollection fec_l2(order - 1, dim);
|
||||
ParFiniteElementSpace fespace_nd(&pmesh, &fec_nd);
|
||||
ParFiniteElementSpace fespace_rt(&pmesh, &fec_rt);
|
||||
ParFiniteElementSpace fespace_ndp(&pmesh, &fec_ndp);
|
||||
ParFiniteElementSpace fespace_rtp(&pmesh, &fec_rtp);
|
||||
ParFiniteElementSpace fespace_l2(&pmesh, &fec_l2);
|
||||
HYPRE_Int size = fespace_nd.GlobalTrueVSize() +
|
||||
fespace_ndp.GlobalTrueVSize();
|
||||
if (myid == 0)
|
||||
{
|
||||
cout << "Number of unknowns: " << size << endl;
|
||||
}
|
||||
|
||||
// 7. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
|
||||
// element space. The first corresponds to the curl curl, while the second
|
||||
// is a simple mass matrix needed on the right hand side of the
|
||||
// generalized eigenvalue problem below. The boundary conditions are
|
||||
// implemented by marking all the boundary attributes from the mesh as
|
||||
// essential. The corresponding degrees of freedom are eliminated with
|
||||
// special values on the diagonal to shift the Dirichlet eigenvalues out
|
||||
// of the computational range. After serial and parallel assembly we
|
||||
// extract the corresponding parallel matrices A and M.
|
||||
FunctionCoefficient rhoCoef(rhoFunc);
|
||||
FunctionCoefficient rhoInvCoef(rhoInvFunc);
|
||||
Array<int> ess_bdr_rz(pmesh.bdr_attributes.Size());
|
||||
Array<int> ess_bdr_phi(pmesh.bdr_attributes.Size());
|
||||
ess_bdr_rz = 1; ess_bdr_rz[3] = 0;
|
||||
ess_bdr_phi = 1;
|
||||
|
||||
ParBilinearForm *arz = new ParBilinearForm(&fespace_nd);
|
||||
arz->AddDomainIntegrator(new CurlCurlIntegrator(rhoCoef));
|
||||
arz->Assemble();
|
||||
arz->EliminateEssentialBCDiag(ess_bdr_rz, 1.0);
|
||||
arz->Finalize();
|
||||
|
||||
ParBilinearForm *mrz = new ParBilinearForm(&fespace_nd);
|
||||
mrz->AddDomainIntegrator(new VectorFEMassIntegrator(rhoCoef));
|
||||
mrz->Assemble();
|
||||
// shift the eigenvalue corresponding to eliminated dofs to a large value
|
||||
mrz->EliminateEssentialBCDiag(ess_bdr_rz, numeric_limits<double>::min());
|
||||
mrz->Finalize();
|
||||
|
||||
ParBilinearForm *aphi = new ParBilinearForm(&fespace_ndp);
|
||||
aphi->AddDomainIntegrator(new MassIntegrator(rhoInvCoef));
|
||||
aphi->AddDomainIntegrator(new DiffusionIntegrator(rhoCoef));
|
||||
aphi->Assemble();
|
||||
aphi->EliminateEssentialBCDiag(ess_bdr_phi, 1.0);
|
||||
aphi->Finalize();
|
||||
|
||||
ParBilinearForm *mphi = new ParBilinearForm(&fespace_ndp);
|
||||
mphi->AddDomainIntegrator(new MassIntegrator(rhoCoef));
|
||||
mphi->Assemble();
|
||||
// shift the eigenvalue corresponding to eliminated dofs to a large value
|
||||
mphi->EliminateEssentialBCDiag(ess_bdr_phi, numeric_limits<double>::min());
|
||||
mphi->Finalize();
|
||||
|
||||
HypreParMatrix *Arz = arz->ParallelAssemble();
|
||||
HypreParMatrix *Mrz = mrz->ParallelAssemble();
|
||||
|
||||
HypreParMatrix *Aphi = aphi->ParallelAssemble();
|
||||
HypreParMatrix *Mphi = mphi->ParallelAssemble();
|
||||
|
||||
delete arz;
|
||||
delete mrz;
|
||||
delete aphi;
|
||||
delete mphi;
|
||||
|
||||
// 8. Define and configure the AME eigensolver and the AMS preconditioner for
|
||||
// A to be used within the solver. Set the matrices which define the
|
||||
// generalized eigenproblem A x = lambda M x.
|
||||
HypreAMS *ams = new HypreAMS(*Arz,&fespace_nd);
|
||||
ams->SetPrintLevel(0);
|
||||
ams->SetSingularProblem();
|
||||
|
||||
HypreAME *ame = new HypreAME(MPI_COMM_WORLD);
|
||||
ame->SetNumModes(nev);
|
||||
ame->SetPreconditioner(*ams);
|
||||
ame->SetMaxIter(100);
|
||||
ame->SetTol(1e-8);
|
||||
ame->SetPrintLevel(1);
|
||||
ame->SetMassMatrix(*Mrz);
|
||||
ame->SetOperator(*Arz);
|
||||
|
||||
HypreBoomerAMG *amg = new HypreBoomerAMG(*Aphi);
|
||||
amg->SetPrintLevel(0);
|
||||
|
||||
HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
|
||||
lobpcg->SetNumModes(nev);
|
||||
// lobpcg->SetRandomSeed(seed);
|
||||
lobpcg->SetPreconditioner(*amg);
|
||||
lobpcg->SetMaxIter(200);
|
||||
lobpcg->SetTol(1e-8);
|
||||
lobpcg->SetPrecondUsageMode(1);
|
||||
lobpcg->SetPrintLevel(1);
|
||||
lobpcg->SetMassMatrix(*Mphi);
|
||||
lobpcg->SetOperator(*Aphi);
|
||||
|
||||
// 9. Compute the eigenmodes and extract the array of eigenvalues. Define a
|
||||
// parallel grid function to represent each of the eigenmodes returned by
|
||||
// the solver. Also, define a discrete curl operator.
|
||||
Array<double> eigenvalues_rz;
|
||||
Array<double> eigenvalues_phi;
|
||||
|
||||
ame->Solve();
|
||||
ame->GetEigenvalues(eigenvalues_rz);
|
||||
|
||||
lobpcg->Solve();
|
||||
lobpcg->GetEigenvalues(eigenvalues_phi);
|
||||
|
||||
Array<double> eigenvalues;
|
||||
|
||||
ParGridFunction x_rz(&fespace_nd);
|
||||
ParGridFunction x_phi(&fespace_ndp);
|
||||
ParGridFunction dx_rz(&fespace_rt);
|
||||
ParGridFunction dx_phi(&fespace_rtp);
|
||||
|
||||
ParDiscreteLinearOperator curl(&fespace_nd, &fespace_rtp);
|
||||
curl.AddDomainInterpolator(new CurlInterpolator());
|
||||
curl.Assemble();
|
||||
curl.Finalize();
|
||||
|
||||
ParDiscreteLinearOperator curl2(&fespace_ndp, &fespace_rt);
|
||||
curl2.AddDomainInterpolator(new CurlInterpolator());
|
||||
curl2.Assemble();
|
||||
curl2.Finalize();
|
||||
|
||||
// This is one workaround for GLVis limitations
|
||||
ParGridFunction dx_l2(&fespace_l2);
|
||||
|
||||
GridFunctionCoefficient dxCoef(&dx_phi);
|
||||
|
||||
if ( myid == 0 )
|
||||
{
|
||||
cout << "\nRelative error in eigenvalues of RZ modes:\n";
|
||||
for (int i=0; i<nev; i++)
|
||||
{
|
||||
double lambda =
|
||||
pow(J0z[mode_nums_rz[2*i]], 2) +
|
||||
pow(M_PI * mode_nums_rz[2*i+1], 2);
|
||||
cout << "Lambda " << i+1 << '/' << nev << " = " << eigenvalues_rz[i]
|
||||
<< ", rel err = " << fabs(eigenvalues_rz[i] - lambda) / lambda
|
||||
<< endl;
|
||||
}
|
||||
cout << endl;
|
||||
cout << "\nRelative error in eigenvalues of Phi modes:\n";
|
||||
for (int i=0; i<nev; i++)
|
||||
{
|
||||
double lambda =
|
||||
pow(J1z[mode_nums_phi[2*i]], 2) +
|
||||
pow(M_PI * mode_nums_phi[2*i+1], 2);
|
||||
cout << "Lambda " << i+1 << '/' << nev << " = " << eigenvalues_phi[i]
|
||||
<< ", rel err = " << fabs(eigenvalues_phi[i] - lambda) / lambda
|
||||
<< endl;
|
||||
}
|
||||
cout << endl;
|
||||
}
|
||||
|
||||
// 10. Save the refined mesh and the modes in parallel. This output can be
|
||||
// viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
|
||||
{
|
||||
ostringstream mesh_name, mode_name;
|
||||
mesh_name << "mesh." << setfill('0') << setw(6) << myid;
|
||||
|
||||
ofstream mesh_ofs(mesh_name.str().c_str());
|
||||
mesh_ofs.precision(8);
|
||||
pmesh.Print(mesh_ofs);
|
||||
|
||||
for (int i=0; i<nev; i++)
|
||||
{
|
||||
// convert eigenvector from HypreParVector to ParGridFunction
|
||||
x_rz = ame->GetEigenvector(i);
|
||||
|
||||
mode_name << "mode_rz_" << setfill('0') << setw(2) << i << "."
|
||||
<< setfill('0') << setw(6) << myid;
|
||||
|
||||
ofstream mode_rz_ofs(mode_name.str().c_str());
|
||||
mode_rz_ofs.precision(8);
|
||||
x_rz.Save(mode_rz_ofs);
|
||||
mode_name.str("");
|
||||
|
||||
// convert eigenvector from HypreParVector to ParGridFunction
|
||||
x_phi = lobpcg->GetEigenvector(i);
|
||||
|
||||
mode_name << "mode_phi_" << setfill('0') << setw(2) << i << "."
|
||||
<< setfill('0') << setw(6) << myid;
|
||||
|
||||
ofstream mode_phi_ofs(mode_name.str().c_str());
|
||||
mode_phi_ofs.precision(8);
|
||||
x_phi.Save(mode_phi_ofs);
|
||||
mode_name.str("");
|
||||
}
|
||||
}
|
||||
|
||||
// 11. Send the solution by socket to a GLVis server.
|
||||
if (visualization)
|
||||
{
|
||||
char vishost[] = "localhost";
|
||||
int visport = 19916;
|
||||
socketstream mode_rz_sock(vishost, visport);
|
||||
socketstream curl_rz_sock(vishost, visport);
|
||||
socketstream mode_phi_sock(vishost, visport);
|
||||
socketstream curl_phi_sock(vishost, visport);
|
||||
mode_rz_sock.precision(8);
|
||||
curl_rz_sock.precision(8);
|
||||
mode_phi_sock.precision(8);
|
||||
curl_phi_sock.precision(8);
|
||||
|
||||
int irz = 0;
|
||||
int iphi = 0;
|
||||
|
||||
for (int i=0; i<nev; i++)
|
||||
{
|
||||
if (mode_type[i] == 0)
|
||||
{
|
||||
if ( myid == 0 )
|
||||
{
|
||||
cout << "Eigenmode " << i+1 << '/' << nev
|
||||
<< ", Lambda = " << eigenvalues_rz[irz] << endl;
|
||||
}
|
||||
|
||||
// convert eigenvector from HypreParVector to ParGridFunction
|
||||
x_rz = ame->GetEigenvector(irz);
|
||||
|
||||
curl.Mult(x_rz, dx_phi);
|
||||
dx_l2.ProjectCoefficient(dxCoef);
|
||||
|
||||
mode_rz_sock << "parallel " << num_procs << " " << myid << "\n"
|
||||
<< "solution\n" << pmesh << x_rz << flush
|
||||
<< "window_title 'Eigenmode " << i+1 << '/' << nev
|
||||
<< ", Lambda = " << eigenvalues_rz[irz] << "'\n";
|
||||
if (irz == 0)
|
||||
{
|
||||
mode_rz_sock << "keys vvv\n";
|
||||
}
|
||||
mode_rz_sock << flush;
|
||||
|
||||
curl_rz_sock << "parallel " << num_procs << " " << myid << "\n"
|
||||
<< "solution\n" << pmesh << dx_l2 << flush
|
||||
<< "window_title 'Curl of Eigenmode " << i+1
|
||||
<< '/' << nev
|
||||
<< ", Lambda = " << eigenvalues_rz[irz] << "' "
|
||||
<< "window_geometry 400 0 400 350\n" << flush;
|
||||
|
||||
irz++;
|
||||
}
|
||||
else
|
||||
{
|
||||
if ( myid == 0 )
|
||||
{
|
||||
cout << "Eigenmode " << i+1 << '/' << nev
|
||||
<< ", Lambda = " << eigenvalues_phi[iphi] << endl;
|
||||
}
|
||||
|
||||
// convert eigenvector from HypreParVector to ParGridFunction
|
||||
x_phi = lobpcg->GetEigenvector(iphi);
|
||||
|
||||
curl2.Mult(x_phi, dx_rz);
|
||||
|
||||
mode_phi_sock << "parallel " << num_procs << " " << myid << "\n"
|
||||
<< "solution\n" << pmesh << x_phi << flush
|
||||
<< "window_title 'Eigenmode " << i+1 << '/' << nev
|
||||
<< ", Lambda = " << eigenvalues_phi[iphi] << "' "
|
||||
<< "window_geometry 0 375 400 350\n"
|
||||
<< flush;
|
||||
|
||||
curl_phi_sock << "parallel " << num_procs << " " << myid << "\n"
|
||||
<< "solution\n" << pmesh << dx_rz << flush
|
||||
<< "window_title 'Curl of Eigenmode " << i+1
|
||||
<< '/' << nev
|
||||
<< ", Lambda = " << eigenvalues_phi[iphi] << "' "
|
||||
<< "window_geometry 400 375 400 350\n";
|
||||
if (iphi == 0)
|
||||
{
|
||||
curl_phi_sock << "keys vvv\n";
|
||||
}
|
||||
curl_phi_sock << flush;
|
||||
|
||||
iphi++;
|
||||
}
|
||||
char c;
|
||||
if (myid == 0)
|
||||
{
|
||||
cout << "press (q)uit or (c)ontinue --> " << flush;
|
||||
cin >> c;
|
||||
}
|
||||
MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
|
||||
|
||||
if (c != 'c')
|
||||
{
|
||||
break;
|
||||
}
|
||||
}
|
||||
mode_rz_sock.close();
|
||||
curl_rz_sock.close();
|
||||
mode_phi_sock.close();
|
||||
curl_phi_sock.close();
|
||||
}
|
||||
/*
|
||||
if (visualization)
|
||||
{
|
||||
char vishost[] = "localhost";
|
||||
int visport = 19916;
|
||||
socketstream mode_sock(vishost, visport);
|
||||
mode_sock.precision(8);
|
||||
socketstream curl_sock(vishost, visport);
|
||||
curl_sock.precision(8);
|
||||
|
||||
for (int i=0; i<nev; i++)
|
||||
{
|
||||
if ( myid == 0 )
|
||||
{
|
||||
cout << "Eigenmode " << i+1 << '/' << nev
|
||||
<< ", Lambda = " << eigenvalues_phi[i] << endl;
|
||||
}
|
||||
|
||||
// convert eigenvector from HypreParVector to ParGridFunction
|
||||
x_phi = lobpcg->GetEigenvector(i);
|
||||
|
||||
curl2.Mult(x_phi, dx_rz);
|
||||
// dx_l2.ProjectCoefficient(dxCoef);
|
||||
|
||||
mode_sock << "parallel " << num_procs << " " << myid << "\n"
|
||||
<< "solution\n" << pmesh << x_phi << flush
|
||||
<< "window_title 'Eigenmode " << i+1 << '/' << nev
|
||||
<< ", Lambda = " << eigenvalues_phi[i] << "' "
|
||||
<< "keys vvv\n" << flush;
|
||||
|
||||
// Limitations in the GridFunction and GLVis prevent this from working
|
||||
curl_sock << "parallel " << num_procs << " " << myid << "\n"
|
||||
<< "solution\n" << pmesh << dx_rz << flush
|
||||
<< "window_title 'Curl of Eigenmode " << i+1 << '/' << nev
|
||||
<< ", Lambda = " << eigenvalues_phi[i] << "' "
|
||||
<< "window_geometry 400 0 400 350\n" << flush;
|
||||
|
||||
char c;
|
||||
if (myid == 0)
|
||||
{
|
||||
cout << "press (q)uit or (c)ontinue --> " << flush;
|
||||
cin >> c;
|
||||
}
|
||||
MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
|
||||
|
||||
if (c != 'c')
|
||||
{
|
||||
break;
|
||||
}
|
||||
}
|
||||
mode_sock.close();
|
||||
}
|
||||
*/
|
||||
// 12. Free the used memory.
|
||||
delete ame;
|
||||
delete ams;
|
||||
delete lobpcg;
|
||||
delete amg;
|
||||
delete Mrz;
|
||||
delete Arz;
|
||||
delete Mphi;
|
||||
delete Aphi;
|
||||
|
||||
MPI_Finalize();
|
||||
|
||||
return 0;
|
||||
}
|
||||
@@ -1,336 +0,0 @@
|
||||
// MFEM Example 13-cyl - Parallel Version
|
||||
//
|
||||
// Compile with: make ex13p-cyl
|
||||
//
|
||||
// Sample runs: mpirun -np 4 ex13p-cyl
|
||||
// mpirun -np 4 ex13p-cyl -o 2
|
||||
// mpirun -np 4 ex13p-cyl -o 2 -e 0
|
||||
//
|
||||
// Description: This example code demonstrates the use of MFEM to solve the
|
||||
// Maxwell (electromagnetic) eigenvalue problem on an
|
||||
// axisymmetric domain. The eigenvalue problem:
|
||||
// curl curl E = lambda E
|
||||
// with homogeneous Dirichlet boundary conditions E x n = 0 is
|
||||
// solved on a cylindrical domain by meshing only a rectangle in
|
||||
// the rho, z plane. In cylindrical coordinates the weak form of
|
||||
// the eigenvalue problem is given by:
|
||||
// (rho Curl(u), Curl(v)) = lambda (rho u, v)
|
||||
//
|
||||
// We compute the five lowest nonzero eigenmodes by discretizing
|
||||
// the curl curl operator using a Nedelec FE space of the
|
||||
// specified order and compare to the known values. Because the
|
||||
// eigenvalue spectrum of a domain is unique this provides a
|
||||
// reliable test that the axisymmetric domain is being faithfully
|
||||
// characterized.
|
||||
//
|
||||
// The example highlights the use of specialized coefficients
|
||||
// with existing operators to mimic axisymmetric domains. The
|
||||
// curl of each eigenmode is also computed and displayed to
|
||||
// ilustrate that no special steps need to be taken to compute
|
||||
// the curl in this coordinate system.
|
||||
//
|
||||
// We recommend viewing examples 13 and 11-cyl before viewing this
|
||||
// example.
|
||||
|
||||
#include "mfem.hpp"
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
|
||||
using namespace std;
|
||||
using namespace mfem;
|
||||
|
||||
// Zeros of Bessel function J_0
|
||||
static double J0z[] = {2.40482555769577,
|
||||
5.52007811028631,
|
||||
8.65372791291101,
|
||||
11.7915344390143
|
||||
};
|
||||
|
||||
// Modes numbers in the rho and z directions for the first five eigenmodes
|
||||
static int mode_nums[] = {0, 0,
|
||||
0, 1,
|
||||
1, 0,
|
||||
1, 1,
|
||||
0, 2
|
||||
};
|
||||
|
||||
double rhoFunc(const Vector &x)
|
||||
{
|
||||
return x[0];
|
||||
}
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
// 1. Initialize MPI.
|
||||
int num_procs, myid;
|
||||
MPI_Init(&argc, &argv);
|
||||
MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
|
||||
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
|
||||
|
||||
// 2. Parse command-line options.
|
||||
int dim = 2;
|
||||
int nr = 1;
|
||||
int nz = 1;
|
||||
int el_type_flag = 1;
|
||||
Element::Type el_type;
|
||||
int ser_ref_levels = 2;
|
||||
int par_ref_levels = 1;
|
||||
int order = 1;
|
||||
int nev = 5;
|
||||
bool visualization = 1;
|
||||
|
||||
OptionsParser args(argc, argv);
|
||||
args.AddOption(&nz, "-nz", "--num-elements-z",
|
||||
"Number of elements in z-direction.");
|
||||
args.AddOption(&nr, "-nr", "--num-elements-rho",
|
||||
"Number of elements in radial direction.");
|
||||
args.AddOption(&el_type_flag, "-e", "--element-type",
|
||||
"Element type: 0 - Triangle, 1 - Quadrilateral.");
|
||||
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",
|
||||
"Finite element order (polynomial degree) or -1 for"
|
||||
" isoparametric space.");
|
||||
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
|
||||
"--no-visualization",
|
||||
"Enable or disable GLVis visualization.");
|
||||
args.Parse();
|
||||
if (!args.Good())
|
||||
{
|
||||
if (myid == 0)
|
||||
{
|
||||
args.PrintUsage(cout);
|
||||
}
|
||||
MPI_Finalize();
|
||||
return 1;
|
||||
}
|
||||
if (myid == 0)
|
||||
{
|
||||
args.PrintOptions(cout);
|
||||
}
|
||||
|
||||
// The output mesh could be quadrilaterals or triangles
|
||||
el_type = (el_type_flag == 0) ? Element::TRIANGLE : Element::QUADRILATERAL;
|
||||
if (el_type != Element::TRIANGLE && el_type != Element::QUADRILATERAL)
|
||||
{
|
||||
cout << "Unsupported element type" << endl;
|
||||
exit(1);
|
||||
}
|
||||
|
||||
// 3. Prepare a rectangular mesh with the desired dimensions and element
|
||||
// type. Other 2D meshes could be used but then we couldn't check the
|
||||
// eigenvalues.
|
||||
ParMesh pmesh;
|
||||
{
|
||||
Mesh mesh = Mesh::MakeCartesian2D(nr, nz, el_type);
|
||||
|
||||
// 4. Refine the serial mesh on all processors to increase the resolution. In
|
||||
// this example we do 'ref_levels' of uniform refinement (2 by default, or
|
||||
// specified on the command line with -rs).
|
||||
for (int lev = 0; lev < ser_ref_levels; lev++)
|
||||
{
|
||||
mesh.UniformRefinement();
|
||||
}
|
||||
|
||||
// 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
|
||||
// this mesh further in parallel to increase the resolution (1 time by
|
||||
// default, or specified on the command line with -rp). Once the parallel
|
||||
// mesh is defined, the serial mesh can be deleted.
|
||||
pmesh = ParMesh(MPI_COMM_WORLD, mesh);
|
||||
for (int lev = 0; lev < par_ref_levels; lev++)
|
||||
{
|
||||
pmesh.UniformRefinement();
|
||||
}
|
||||
pmesh.ReorientTetMesh();
|
||||
}
|
||||
|
||||
// 6. Define a parallel finite element space on the parallel mesh. Here we
|
||||
// use the Nedelec finite elements (ND) of the specified order. We also
|
||||
// create an L2 space to represent the z-component of the curl of the
|
||||
// modes.
|
||||
ND_FECollection fec_nd(order, dim);
|
||||
L2_FECollection fec_rt(order - 1, dim,
|
||||
BasisType::GaussLegendre, FiniteElement::INTEGRAL);
|
||||
L2_FECollection fec_l2(order - 1, dim);
|
||||
ParFiniteElementSpace fespace_nd(&pmesh, &fec_nd);
|
||||
ParFiniteElementSpace fespace_rt(&pmesh, &fec_rt);
|
||||
ParFiniteElementSpace fespace_l2(&pmesh, &fec_l2);
|
||||
HYPRE_Int size = fespace_nd.GlobalTrueVSize();
|
||||
if (myid == 0)
|
||||
{
|
||||
cout << "Number of unknowns: " << size << endl;
|
||||
}
|
||||
|
||||
// 7. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
|
||||
// element space. The first corresponds to the curl curl, while the second
|
||||
// is a simple mass matrix needed on the right hand side of the
|
||||
// generalized eigenvalue problem below. The boundary conditions are
|
||||
// implemented by marking all the boundary attributes from the mesh as
|
||||
// essential. The corresponding degrees of freedom are eliminated with
|
||||
// special values on the diagonal to shift the Dirichlet eigenvalues out
|
||||
// of the computational range. After serial and parallel assembly we
|
||||
// extract the corresponding parallel matrices A and M.
|
||||
FunctionCoefficient rhoCoef(rhoFunc);
|
||||
Array<int> ess_bdr(pmesh.bdr_attributes.Size());
|
||||
ess_bdr = 1;
|
||||
ess_bdr[3] = 0;
|
||||
|
||||
ParBilinearForm *a = new ParBilinearForm(&fespace_nd);
|
||||
a->AddDomainIntegrator(new CurlCurlIntegrator(rhoCoef));
|
||||
a->Assemble();
|
||||
a->EliminateEssentialBCDiag(ess_bdr, 1.0);
|
||||
a->Finalize();
|
||||
|
||||
ParBilinearForm *m = new ParBilinearForm(&fespace_nd);
|
||||
m->AddDomainIntegrator(new VectorFEMassIntegrator(rhoCoef));
|
||||
m->Assemble();
|
||||
// shift the eigenvalue corresponding to eliminated dofs to a large value
|
||||
m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
|
||||
m->Finalize();
|
||||
|
||||
HypreParMatrix *A = a->ParallelAssemble();
|
||||
HypreParMatrix *M = m->ParallelAssemble();
|
||||
|
||||
delete a;
|
||||
delete m;
|
||||
|
||||
// 8. Define and configure the AME eigensolver and the AMS preconditioner for
|
||||
// A to be used within the solver. Set the matrices which define the
|
||||
// generalized eigenproblem A x = lambda M x.
|
||||
HypreAMS *ams = new HypreAMS(*A,&fespace_nd);
|
||||
ams->SetPrintLevel(0);
|
||||
ams->SetSingularProblem();
|
||||
|
||||
HypreAME *ame = new HypreAME(MPI_COMM_WORLD);
|
||||
ame->SetNumModes(nev);
|
||||
ame->SetPreconditioner(*ams);
|
||||
ame->SetMaxIter(100);
|
||||
ame->SetTol(1e-8);
|
||||
ame->SetPrintLevel(1);
|
||||
ame->SetMassMatrix(*M);
|
||||
ame->SetOperator(*A);
|
||||
|
||||
// 9. Compute the eigenmodes and extract the array of eigenvalues. Define a
|
||||
// parallel grid function to represent each of the eigenmodes returned by
|
||||
// the solver. Also, define a discrete curl operator.
|
||||
Array<double> eigenvalues;
|
||||
ame->Solve();
|
||||
ame->GetEigenvalues(eigenvalues);
|
||||
ParGridFunction x(&fespace_nd);
|
||||
ParGridFunction dx(&fespace_rt);
|
||||
|
||||
ParDiscreteLinearOperator curl(&fespace_nd, &fespace_rt);
|
||||
curl.AddDomainInterpolator(new CurlInterpolator());
|
||||
curl.Assemble();
|
||||
curl.Finalize();
|
||||
|
||||
// This is one workaround for GLVis limitations
|
||||
ParGridFunction dx_l2(&fespace_l2);
|
||||
|
||||
GridFunctionCoefficient dxCoef(&dx);
|
||||
|
||||
if ( myid == 0 )
|
||||
{
|
||||
cout << "\nRelative error in eigenvalues:\n";
|
||||
for (int i=0; i<nev; i++)
|
||||
{
|
||||
double lambda =
|
||||
pow(J0z[mode_nums[2*i]], 2) +
|
||||
pow(M_PI * mode_nums[2*i+1], 2);
|
||||
cout << "Lambda " << i+1 << '/' << nev << " = " << eigenvalues[i]
|
||||
<< ", rel err = " << fabs(eigenvalues[i] - lambda) / lambda
|
||||
<< endl;
|
||||
}
|
||||
cout << endl;
|
||||
}
|
||||
|
||||
// 10. Save the refined mesh and the modes in parallel. This output can be
|
||||
// viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
|
||||
{
|
||||
ostringstream mesh_name, mode_name;
|
||||
mesh_name << "mesh." << setfill('0') << setw(6) << myid;
|
||||
|
||||
ofstream mesh_ofs(mesh_name.str().c_str());
|
||||
mesh_ofs.precision(8);
|
||||
pmesh.Print(mesh_ofs);
|
||||
|
||||
for (int i=0; i<nev; i++)
|
||||
{
|
||||
// convert eigenvector from HypreParVector to ParGridFunction
|
||||
x = ame->GetEigenvector(i);
|
||||
|
||||
mode_name << "mode_" << setfill('0') << setw(2) << i << "."
|
||||
<< setfill('0') << setw(6) << myid;
|
||||
|
||||
ofstream mode_ofs(mode_name.str().c_str());
|
||||
mode_ofs.precision(8);
|
||||
x.Save(mode_ofs);
|
||||
mode_name.str("");
|
||||
}
|
||||
}
|
||||
|
||||
// 11. Send the solution by socket to a GLVis server.
|
||||
if (visualization)
|
||||
{
|
||||
char vishost[] = "localhost";
|
||||
int visport = 19916;
|
||||
socketstream mode_sock(vishost, visport);
|
||||
mode_sock.precision(8);
|
||||
socketstream curl_sock(vishost, visport);
|
||||
curl_sock.precision(8);
|
||||
|
||||
for (int i=0; i<nev; i++)
|
||||
{
|
||||
if ( myid == 0 )
|
||||
{
|
||||
cout << "Eigenmode " << i+1 << '/' << nev
|
||||
<< ", Lambda = " << eigenvalues[i] << endl;
|
||||
}
|
||||
|
||||
// convert eigenvector from HypreParVector to ParGridFunction
|
||||
x = ame->GetEigenvector(i);
|
||||
|
||||
curl.Mult(x, dx);
|
||||
dx_l2.ProjectCoefficient(dxCoef);
|
||||
|
||||
mode_sock << "parallel " << num_procs << " " << myid << "\n"
|
||||
<< "solution\n" << pmesh << x << flush
|
||||
<< "window_title 'Eigenmode " << i+1 << '/' << nev
|
||||
<< ", Lambda = " << eigenvalues[i] << "' "
|
||||
<< "keys vvv\n" << flush;
|
||||
|
||||
// Limitations in the GridFunction and GLVis prevent this from working
|
||||
curl_sock << "parallel " << num_procs << " " << myid << "\n"
|
||||
<< "solution\n" << pmesh << dx_l2 << flush
|
||||
<< "window_title 'Curl of Eigenmode " << i+1 << '/' << nev
|
||||
<< ", Lambda = " << eigenvalues[i] << "' "
|
||||
<< "window_geometry 400 0 400 350\n" << flush;
|
||||
|
||||
char c;
|
||||
if (myid == 0)
|
||||
{
|
||||
cout << "press (q)uit or (c)ontinue --> " << flush;
|
||||
cin >> c;
|
||||
}
|
||||
MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
|
||||
|
||||
if (c != 'c')
|
||||
{
|
||||
break;
|
||||
}
|
||||
}
|
||||
mode_sock.close();
|
||||
}
|
||||
|
||||
// 12. Free the used memory.
|
||||
delete ame;
|
||||
delete ams;
|
||||
delete M;
|
||||
delete A;
|
||||
|
||||
MPI_Finalize();
|
||||
|
||||
return 0;
|
||||
}
|
||||
@@ -1,975 +0,0 @@
|
||||
// MFEM Example 1 Ortho - Parallel Version
|
||||
//
|
||||
// Compile with: make ex1p
|
||||
//
|
||||
// Sample runs: mpirun -np 4 ex1p-orth
|
||||
// mpirun -np 4 ex1p-orth -c 2
|
||||
// mpirun -np 4 ex1p-orth -c 3
|
||||
// mpirun -np 4 ex1p-orth -c 4 -rs 4
|
||||
// mpirun -np 4 ex1p-orth -c 6
|
||||
// mpirun -np 4 ex1p-orth -c 7
|
||||
// mpirun -np 4 ex1p-orth -c 8
|
||||
// mpirun -np 4 ex1p-orth -c 9
|
||||
// mpirun -np 4 ex1p-orth -c 10 -rs 4
|
||||
// mpirun -np 4 ex1p-orth -c 11 -n2 2 -rs 3
|
||||
//
|
||||
// Description: This example code demonstrates the use of MFEM to define a
|
||||
// simple finite element discretization of the Laplace problem
|
||||
// -Delta u = 1 with homogeneous Dirichlet boundary conditions in
|
||||
// a variety of orthogonal coordinate systems. The discretization
|
||||
// is identical to that used in example 1 but here we use
|
||||
// non-trivial coefficients in the Laplace operator and the right-
|
||||
// hand-side vector to mimic a curvilinear coordinate system. We
|
||||
// also transform the mesh and solve the standard Laplace problem
|
||||
// on the transformed mesh to compare the solutions.
|
||||
//
|
||||
// The example highlights the use of standard differential
|
||||
// operators to mimic the behavior of more exotic operators
|
||||
// derived from coordinate transformations.
|
||||
//
|
||||
// We recommend viewing Example 1 and Example 11-cyl before
|
||||
// viewing this example.
|
||||
//
|
||||
// Note: the notation used in this code comes from the Wikipedia
|
||||
// page https://en.wikipedia.org/wiki/Orthogonal_coordinates.
|
||||
// There are, however, minor differences made to ensure that
|
||||
// we use right-handed coordinate systems in all cases.
|
||||
|
||||
#include "mfem.hpp"
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
|
||||
using namespace std;
|
||||
using namespace mfem;
|
||||
|
||||
// Enumeration listing all supported 2D orthogonal coordinate systems
|
||||
enum CoordSys {POLAR = 1, PARABOLIC_CYL, ELLIPTIC, BIPOLAR,
|
||||
CYLINDRICAL, SPHERICAL, PARABOLIC, PROLATE_SPHEROIDAL,
|
||||
OBLATE_SPHEROIDAL, TOROIDAL, BISPHERICAL
|
||||
};
|
||||
|
||||
static CoordSys coords_ = (CoordSys)1;
|
||||
static double q1_min_ = NAN;
|
||||
static double q1_max_ = NAN;
|
||||
static double q2_min_ = NAN;
|
||||
static double q2_max_ = NAN;
|
||||
static double a_ = 1.0;
|
||||
|
||||
// Set default values for coordinate ranges q1_min_, q1_max_, q2_min_,
|
||||
// and q2_max_ based on the selected coordinate system, coords_.
|
||||
void SetRanges();
|
||||
|
||||
// Shift the mesh so that the origin is at (q1_min_, q2_min_)
|
||||
void trans1(const Vector &u, Vector &x)
|
||||
{
|
||||
x.SetSize(2);
|
||||
x[0] = u[0] + q1_min_;
|
||||
x[1] = u[1] + q2_min_;
|
||||
}
|
||||
|
||||
// Apply conformal mapping from cartesian coordinates to the
|
||||
// orthogonal coordinate system specified by coords_.
|
||||
void trans(const Vector &u, Vector &x);
|
||||
|
||||
// Returns one of the three coordinate scale factors h_i describing
|
||||
// the orthogonal coordinate system.
|
||||
class OrthoCoef : public Coefficient
|
||||
{
|
||||
private:
|
||||
int ind_;
|
||||
|
||||
public:
|
||||
OrthoCoef(int index) : ind_(index) {}
|
||||
|
||||
virtual double Eval(ElementTransformation &T,
|
||||
const IntegrationPoint &ip);
|
||||
};
|
||||
|
||||
// Integration weight coefficient h_1 * h_2 * h_3
|
||||
class OrthoWeightCoef : public Coefficient
|
||||
{
|
||||
private:
|
||||
Coefficient &h1Coef_;
|
||||
Coefficient &h2Coef_;
|
||||
Coefficient &h3Coef_;
|
||||
|
||||
public:
|
||||
OrthoWeightCoef(Coefficient &h1Coef,
|
||||
Coefficient &h2Coef,
|
||||
Coefficient &h3Coef)
|
||||
: h1Coef_(h1Coef),
|
||||
h2Coef_(h2Coef),
|
||||
h3Coef_(h3Coef) {}
|
||||
|
||||
virtual double Eval(ElementTransformation &T,
|
||||
const IntegrationPoint &ip);
|
||||
};
|
||||
|
||||
// Matrix-valued coefficient appearing in the weak form of the
|
||||
// Laplacian operator.
|
||||
class OrthoMatrixCoef : public MatrixCoefficient
|
||||
{
|
||||
private:
|
||||
Coefficient &h1Coef_;
|
||||
Coefficient &h2Coef_;
|
||||
Coefficient &h3Coef_;
|
||||
|
||||
public:
|
||||
OrthoMatrixCoef(Coefficient &h1Coef,
|
||||
Coefficient &h2Coef,
|
||||
Coefficient &h3Coef)
|
||||
: MatrixCoefficient(2),
|
||||
h1Coef_(h1Coef),
|
||||
h2Coef_(h2Coef),
|
||||
h3Coef_(h3Coef)
|
||||
{}
|
||||
|
||||
virtual void Eval(DenseMatrix &K, ElementTransformation &T,
|
||||
const IntegrationPoint &ip);
|
||||
|
||||
};
|
||||
|
||||
// Radial weight factor to distinguish volumes of revolution from
|
||||
// extruded volumes
|
||||
class RhoCoef : public Coefficient
|
||||
{
|
||||
public:
|
||||
RhoCoef() {}
|
||||
|
||||
virtual double Eval(ElementTransformation &T,
|
||||
const IntegrationPoint &ip)
|
||||
{
|
||||
double u_data[2];
|
||||
Vector u(u_data, 2);
|
||||
T.Transform(ip, u);
|
||||
return u[0];
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
// Matrix-valued radial weight factor to distinguish volumes of
|
||||
// revolution from extruded volumes within the Laplacian operator
|
||||
class RhoMatrixCoef : public MatrixCoefficient
|
||||
{
|
||||
public:
|
||||
RhoMatrixCoef()
|
||||
: MatrixCoefficient(2)
|
||||
{}
|
||||
|
||||
virtual void Eval(DenseMatrix &K, ElementTransformation &T,
|
||||
const IntegrationPoint &ip)
|
||||
{
|
||||
double u_data[2];
|
||||
Vector u(u_data, 2);
|
||||
T.Transform(ip, u);
|
||||
|
||||
K.SetSize(2);
|
||||
K(0,0) = u[0];
|
||||
K(0,1) = 0.0;
|
||||
K(1,0) = 0.0;
|
||||
K(1,1) = u[0];
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
static bool static_cond_ = false;
|
||||
static bool pa_ = false;
|
||||
|
||||
// Ensure that m >= 3 if a periodic mesh has been selected
|
||||
void AdjustDimensions(int &m, int &n, int & rs, int & rp);
|
||||
|
||||
// Setup and solve the Poisson problem with boundary conditions
|
||||
// appropriate to the selected coordinate system.
|
||||
void Poisson(ParMesh &pmesh, ParFiniteElementSpace &fespace,
|
||||
MatrixCoefficient &LCoef, Coefficient &MCoef,
|
||||
ParGridFunction &x);
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
// 1. Initialize MPI.
|
||||
int num_procs, myid;
|
||||
MPI_Init(&argc, &argv);
|
||||
MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
|
||||
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
|
||||
|
||||
// 2. Parse command-line options.
|
||||
int coords = 1;
|
||||
int n1 = 1;
|
||||
int n2 = 1;
|
||||
int el_type_flag = 1;
|
||||
Element::Type el_type;
|
||||
int ser_ref_levels = 2;
|
||||
int par_ref_levels = 1;
|
||||
int morder = 2;
|
||||
int order = 2;
|
||||
const char *device_config = "cpu";
|
||||
bool comp = true;
|
||||
bool discont = false;
|
||||
bool visualization = true;
|
||||
|
||||
OptionsParser args(argc, argv);
|
||||
args.AddOption(&coords, "-c", "--coord-sys",
|
||||
"Coordinate system: 1 - POLAR, 2 - PARABOLIC_CYL, "
|
||||
"3 - ELLIPTIC, 4 - BIPOLAR, 5 - CYLINDRICAL, 6 - SPHERICAL, "
|
||||
"7 - PARABOLIC, 8 - PROLATE_SPHEROIDAL, "
|
||||
"9 - OBLATE_SPHEROIDAL, 10 - TOROIDAL, 11 - BISPHERICAL");
|
||||
args.AddOption(&n1, "-n1", "--num-elements-1",
|
||||
"Number of elements in q1-direction.");
|
||||
args.AddOption(&n2, "-n2", "--num-elements-2",
|
||||
"Number of elements in q2-direction.");
|
||||
args.AddOption(&q1_min_, "-q1-min", "--q1-minimum-1",
|
||||
"Minimum value of q1 coordinate.");
|
||||
args.AddOption(&q1_max_, "-q1-max", "--q1-maximum-1",
|
||||
"Maximum value of q1 coordinate.");
|
||||
args.AddOption(&q2_min_, "-q2-min", "--q2-minimum-1",
|
||||
"Minimum value of q2 coordinate.");
|
||||
args.AddOption(&q2_max_, "-q2-max", "--q2-maximum-1",
|
||||
"Maximum value of q2 coordinate.");
|
||||
args.AddOption(&a_, "-a", "--scale-parameter",
|
||||
"Scale paramter appearing in some of the transformations.");
|
||||
args.AddOption(&el_type_flag, "-e", "--element-type",
|
||||
"Element type: 0 - Triangle, 1 - Quadrilateral.");
|
||||
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(&morder, "-mo", "--mesh-order",
|
||||
"Order (polynomial degree) for the mesh geometry.");
|
||||
args.AddOption(&order, "-o", "--order",
|
||||
"Finite element order (polynomial degree) or -1 for"
|
||||
" isoparametric space.");
|
||||
args.AddOption(&comp, "-comp", "--compare", "-no-comp",
|
||||
"--no-compare", "Compare to standard curved mesh solution.");
|
||||
args.AddOption(&static_cond_, "-sc", "--static-condensation", "-no-sc",
|
||||
"--no-static-condensation", "Enable static condensation.");
|
||||
args.AddOption(&pa_, "-pa", "--partial-assembly", "-no-pa",
|
||||
"--no-partial-assembly", "Enable Partial Assembly.");
|
||||
args.AddOption(&device_config, "-d", "--device",
|
||||
"Device configuration string, see Device::Configure().");
|
||||
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
|
||||
"--no-visualization",
|
||||
"Enable or disable GLVis visualization.");
|
||||
args.Parse();
|
||||
if (!args.Good())
|
||||
{
|
||||
if (myid == 0)
|
||||
{
|
||||
args.PrintUsage(cout);
|
||||
}
|
||||
MPI_Finalize();
|
||||
return 1;
|
||||
}
|
||||
if (myid == 0)
|
||||
{
|
||||
args.PrintOptions(cout);
|
||||
}
|
||||
|
||||
// Cast the user input to the enumerated type and set the appropriate
|
||||
// coordinate ranges.
|
||||
coords_ = (CoordSys)coords;
|
||||
SetRanges();
|
||||
|
||||
// The output mesh could be quadrilaterals or triangles
|
||||
el_type = (el_type_flag == 0) ? Element::TRIANGLE : Element::QUADRILATERAL;
|
||||
if (el_type != Element::TRIANGLE && el_type != Element::QUADRILATERAL)
|
||||
{
|
||||
cout << "Unsupported element type" << endl;
|
||||
exit(1);
|
||||
}
|
||||
|
||||
if (coords_ == BIPOLAR || coords_ == TOROIDAL)
|
||||
{
|
||||
AdjustDimensions(n1, n2, ser_ref_levels, par_ref_levels);
|
||||
}
|
||||
else if (coords_ == POLAR || coords_ == ELLIPTIC)
|
||||
{
|
||||
AdjustDimensions(n2, n1, ser_ref_levels, par_ref_levels);
|
||||
}
|
||||
|
||||
// 3. Enable hardware devices such as GPUs, and programming models such as
|
||||
// CUDA, OCCA, RAJA and OpenMP based on command line options.
|
||||
Device device(device_config);
|
||||
if (myid == 0) { device.Print(); }
|
||||
|
||||
// 3. Prepare a rectangular mesh with the desired dimensions and element
|
||||
// type.
|
||||
Mesh *mesh = new Mesh(n1, n2, el_type, false,
|
||||
q1_max_ - q1_min_, q2_max_ - q2_min_);
|
||||
mesh->Transform(trans1);
|
||||
int dim = mesh->Dimension();
|
||||
|
||||
if (coords_ == POLAR || coords_ == ELLIPTIC || coords_ == BIPOLAR ||
|
||||
coords_ == TOROIDAL)
|
||||
{
|
||||
// 4. Stitch the ends of the mesh together
|
||||
discont = true;
|
||||
mesh->SetCurvature(1, discont, 2, Ordering::byVDIM);
|
||||
Array<int> v2v(mesh->GetNV());
|
||||
for (int i = 0; i < v2v.Size(); i++)
|
||||
{
|
||||
v2v[i] = i;
|
||||
}
|
||||
|
||||
if (coords_ == POLAR || coords_ == ELLIPTIC)
|
||||
{
|
||||
// identify vertices at the extremes of the mesh in the q2 direction
|
||||
for (int i=0; i<n1 + 1; i++)
|
||||
{
|
||||
v2v[v2v.Size() - n1 - 1 + i] = i;
|
||||
}
|
||||
}
|
||||
else if (coords_ == BIPOLAR || coords_ == TOROIDAL)
|
||||
{
|
||||
// identify vertices at the extremes of the mesh in the q1 direction
|
||||
for (int i=0; i<n2 + 1; i++)
|
||||
{
|
||||
v2v[(n1 + 1) * i + n1] = (n1 + 1) * i;
|
||||
}
|
||||
}
|
||||
// renumber elements
|
||||
for (int i = 0; i < mesh->GetNE(); i++)
|
||||
{
|
||||
Element *el = mesh->GetElement(i);
|
||||
int *v = el->GetVertices();
|
||||
int nv = el->GetNVertices();
|
||||
for (int j = 0; j < nv; j++)
|
||||
{
|
||||
v[j] = v2v[v[j]];
|
||||
}
|
||||
}
|
||||
// renumber boundary elements
|
||||
for (int i = 0; i < mesh->GetNBE(); i++)
|
||||
{
|
||||
Element *el = mesh->GetBdrElement(i);
|
||||
int *v = el->GetVertices();
|
||||
int nv = el->GetNVertices();
|
||||
for (int j = 0; j < nv; j++)
|
||||
{
|
||||
v[j] = v2v[v[j]];
|
||||
}
|
||||
}
|
||||
mesh->RemoveUnusedVertices();
|
||||
mesh->RemoveInternalBoundaries();
|
||||
mesh->FinalizeTopology();
|
||||
}
|
||||
|
||||
// 5. Refine the serial mesh on all processors to increase the resolution.
|
||||
{
|
||||
for (int l = 0; l < ser_ref_levels; l++)
|
||||
{
|
||||
mesh->UniformRefinement();
|
||||
}
|
||||
}
|
||||
|
||||
// 6. 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_ortho = new ParMesh(MPI_COMM_WORLD, *mesh);
|
||||
delete mesh;
|
||||
{
|
||||
for (int l = 0; l < par_ref_levels; l++)
|
||||
{
|
||||
pmesh_ortho->UniformRefinement();
|
||||
}
|
||||
}
|
||||
|
||||
// 7. Create a standard curved mesh to describe the same geometry
|
||||
ParMesh *pmesh_curved = new ParMesh(*pmesh_ortho);
|
||||
pmesh_curved->SetCurvature(morder, discont);
|
||||
pmesh_curved->Transform(trans);
|
||||
|
||||
// 8. Define a parallel finite element space on the parallel mesh. Here we
|
||||
// use continuous Lagrange finite elements of the specified order.
|
||||
H1_FECollection fec(order, dim);
|
||||
ParFiniteElementSpace fespace_ortho(pmesh_ortho, &fec);
|
||||
HYPRE_Int size = fespace_ortho.GlobalTrueVSize();
|
||||
if (myid == 0)
|
||||
{
|
||||
cout << "Number of finite element unknowns: " << size << endl;
|
||||
}
|
||||
|
||||
// 9. Declare the coordinate scaling factors and the coefficients
|
||||
// needed to form the mass matrix and Laplacian.
|
||||
OrthoCoef h1Coef(0);
|
||||
OrthoCoef h2Coef(1);
|
||||
OrthoCoef h3Coef(2);
|
||||
OrthoWeightCoef WCoef(h1Coef, h2Coef, h3Coef);
|
||||
OrthoMatrixCoef LCoef(h1Coef, h2Coef, h3Coef);
|
||||
|
||||
// 10. Setup and solve the Poisson problem on the cartesian mesh
|
||||
ParGridFunction x_ortho(&fespace_ortho); x_ortho = 0.0;
|
||||
Poisson(*pmesh_ortho, fespace_ortho, LCoef, WCoef, x_ortho);
|
||||
|
||||
// 11. Save the refined mesh and the solution in parallel. This output can
|
||||
// be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
|
||||
{
|
||||
ostringstream mesh_name, sol_name;
|
||||
mesh_name << "mesh_ortho." << setfill('0') << setw(6) << myid;
|
||||
sol_name << "sol_ortho." << setfill('0') << setw(6) << myid;
|
||||
|
||||
ofstream mesh_ofs(mesh_name.str().c_str());
|
||||
mesh_ofs.precision(8);
|
||||
pmesh_ortho->Print(mesh_ofs);
|
||||
|
||||
ofstream sol_ofs(sol_name.str().c_str());
|
||||
sol_ofs.precision(8);
|
||||
x_ortho.Save(sol_ofs);
|
||||
}
|
||||
|
||||
// 12. Send the solution by socket to a GLVis server.
|
||||
if (visualization)
|
||||
{
|
||||
char vishost[] = "localhost";
|
||||
int visport = 19916;
|
||||
socketstream sol_sock(vishost, visport);
|
||||
sol_sock << "parallel " << num_procs << " " << myid << "\n";
|
||||
sol_sock.precision(8);
|
||||
sol_sock << "solution\n" << *pmesh_ortho << x_ortho << flush
|
||||
<< "window_title 'Straight Mesh'"
|
||||
<< "keys m\n";
|
||||
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
socketstream mix_sol_sock(vishost, visport);
|
||||
mix_sol_sock << "parallel " << num_procs << " " << myid << "\n";
|
||||
mix_sol_sock.precision(8);
|
||||
mix_sol_sock << "solution\n" << *pmesh_curved << x_ortho << flush
|
||||
<< "window_title 'Straight Solution - Curved Mesh' "
|
||||
<< "window_geometry 400 0 400 350"
|
||||
<< "keys m\n";
|
||||
}
|
||||
|
||||
// 13. Compare to a solution computed on the corresponding curved mesh.
|
||||
if (comp)
|
||||
{
|
||||
ParFiniteElementSpace fespace_curved(pmesh_curved, &fec);
|
||||
ParGridFunction x_curved(&fespace_curved); x_curved = 0.0;
|
||||
|
||||
double err = -1.0;
|
||||
GridFunctionCoefficient xCoef(&x_ortho);
|
||||
|
||||
// 14. Setup and solve the Poisson problem on the cartesian mesh
|
||||
if (coords_ == POLAR || coords_ == PARABOLIC_CYL ||
|
||||
coords_ == ELLIPTIC || coords_ == BIPOLAR)
|
||||
{
|
||||
// These coordinate systems can be viewed as truly two-dimensional
|
||||
// or simply extruded into the third dimension and so they require
|
||||
// no special coefficients.
|
||||
DenseMatrix OneMat(2);
|
||||
OneMat = 0.0; OneMat(0,0) = 1.0; OneMat(1,1) = 1.0;
|
||||
MatrixConstantCoefficient OneCoef(OneMat);
|
||||
ConstantCoefficient oneCoef(1.0);
|
||||
Poisson(*pmesh_curved, fespace_curved, OneCoef, oneCoef, x_curved);
|
||||
|
||||
// 15a. Measure the difference in the two solutions using an L2 norm.
|
||||
err = x_curved.ComputeL2Error(xCoef);
|
||||
}
|
||||
else
|
||||
{
|
||||
// The remaining coordinate systems are truly three-dimensional and
|
||||
// involve rotation about the second coordinate axis. Consequently,
|
||||
// they require a radial scale factor both in the mass matrix and the
|
||||
// Laplacian operator.
|
||||
RhoCoef rhoCoef;
|
||||
RhoMatrixCoef RhoCoef;
|
||||
Poisson(*pmesh_curved, fespace_curved, RhoCoef, rhoCoef, x_curved);
|
||||
|
||||
// 15b. Measure the difference in the two solutions using an L2 norm.
|
||||
PowerCoefficient sqrtRhoCoef(rhoCoef, 0.5);
|
||||
ProductCoefficient rxoCoef(sqrtRhoCoef, xCoef);
|
||||
GridFunctionCoefficient xcCoef(&x_curved);
|
||||
ProductCoefficient rxcCoef(sqrtRhoCoef, xcCoef);
|
||||
ParGridFunction rx_curved(&fespace_curved);
|
||||
rx_curved.ProjectCoefficient(rxcCoef);
|
||||
err = rx_curved.ComputeL2Error(rxoCoef);
|
||||
}
|
||||
|
||||
if (myid == 0)
|
||||
{
|
||||
cout << "\n|| u_curved - u_ortho ||_{L^2} = " << err << '\n' << endl;
|
||||
}
|
||||
|
||||
// 15. Save the refined mesh and the solution in parallel. This output can
|
||||
// be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
|
||||
{
|
||||
ostringstream mesh_name, sol_name;
|
||||
mesh_name << "mesh_std." << setfill('0') << setw(6) << myid;
|
||||
sol_name << "sol_std." << setfill('0') << setw(6) << myid;
|
||||
|
||||
ofstream mesh_ofs(mesh_name.str().c_str());
|
||||
mesh_ofs.precision(8);
|
||||
pmesh_curved->Print(mesh_ofs);
|
||||
|
||||
ofstream sol_ofs(sol_name.str().c_str());
|
||||
sol_ofs.precision(8);
|
||||
x_curved.Save(sol_ofs);
|
||||
}
|
||||
|
||||
// 16. Send the solution by socket to a GLVis server.
|
||||
if (visualization)
|
||||
{
|
||||
char vishost[] = "localhost";
|
||||
int visport = 19916;
|
||||
|
||||
socketstream cart_sol_sock(vishost, visport);
|
||||
cart_sol_sock << "parallel " << num_procs << " " << myid << "\n";
|
||||
cart_sol_sock.precision(8);
|
||||
cart_sol_sock << "solution\n" << *pmesh_curved << x_curved << flush
|
||||
<< "window_title 'Curved Mesh' "
|
||||
<< "window_geometry 800 0 400 350"
|
||||
<< "keys m\n";
|
||||
}
|
||||
|
||||
delete pmesh_curved;
|
||||
}
|
||||
|
||||
// 17. Free the used memory.
|
||||
delete pmesh_ortho;
|
||||
|
||||
MPI_Finalize();
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
void AdjustDimensions(int &m, int &n, int & rs, int & rp)
|
||||
{
|
||||
while (m < 3 && rs + rp > 0)
|
||||
{
|
||||
m *= 2;
|
||||
n *= 2;
|
||||
(rs > 0) ? rs-- : rp--;
|
||||
}
|
||||
if (m < 3) { m = 3; }
|
||||
}
|
||||
|
||||
void Poisson(ParMesh &pmesh, ParFiniteElementSpace &fespace,
|
||||
MatrixCoefficient &LCoef, Coefficient &MCoef,
|
||||
ParGridFunction &x)
|
||||
{
|
||||
// 8. Determine the list of true (i.e. parallel conforming) essential
|
||||
// boundary dofs. In this example, the boundary conditions are defined
|
||||
// by marking all the boundary attributes from the mesh as essential
|
||||
// (Dirichlet) and converting them to a list of true dofs.
|
||||
Array<int> ess_tdof_list;
|
||||
if (pmesh.bdr_attributes.Size())
|
||||
{
|
||||
Array<int> ess_bdr(pmesh.bdr_attributes.Max());
|
||||
ess_bdr = 1;
|
||||
if (coords_ == CYLINDRICAL)
|
||||
{
|
||||
ess_bdr[3] = 0;
|
||||
}
|
||||
else if (coords_ == SPHERICAL || coords_ == PROLATE_SPHEROIDAL ||
|
||||
coords_ == OBLATE_SPHEROIDAL)
|
||||
{
|
||||
ess_bdr[0] = 0;
|
||||
ess_bdr[2] = 0;
|
||||
}
|
||||
else if (coords_ == PARABOLIC)
|
||||
{
|
||||
ess_bdr[0] = 0;
|
||||
ess_bdr[3] = 0;
|
||||
}
|
||||
else if (coords_ == BISPHERICAL)
|
||||
{
|
||||
ess_bdr[1] = 0;
|
||||
ess_bdr[3] = 0;
|
||||
}
|
||||
fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
|
||||
}
|
||||
|
||||
// 9. Set up the parallel linear form b(.) which corresponds to the
|
||||
// right-hand side of the FEM linear system, which in this case is
|
||||
// (1,phi_i) where phi_i are the basis functions in fespace.
|
||||
ParLinearForm b(&fespace);
|
||||
b.AddDomainIntegrator(new DomainLFIntegrator(MCoef));
|
||||
b.Assemble();
|
||||
|
||||
// 10. Define the solution vector x as a parallel finite element grid function
|
||||
// corresponding to fespace. Initialize x with initial guess of zero,
|
||||
// which satisfies the boundary conditions.
|
||||
// ParGridFunction x(fespace);
|
||||
// x = 0.0;
|
||||
|
||||
// 11. Set up the parallel bilinear form a(.,.) on the finite element space
|
||||
// corresponding to the Laplacian operator -Delta, by adding the Diffusion
|
||||
// domain integrator.
|
||||
ParBilinearForm a(&fespace);
|
||||
if (pa_) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
|
||||
a.AddDomainIntegrator(new DiffusionIntegrator(LCoef));
|
||||
|
||||
// 12. Assemble the parallel bilinear form and the corresponding linear
|
||||
// system, applying any necessary transformations such as: parallel
|
||||
// assembly, eliminating boundary conditions, applying conforming
|
||||
// constraints for non-conforming AMR, static condensation, etc.
|
||||
if (static_cond_) { a.EnableStaticCondensation(); }
|
||||
a.Assemble();
|
||||
|
||||
OperatorPtr A;
|
||||
Vector B, X;
|
||||
a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
|
||||
|
||||
// 13. Solve the linear system A X = B.
|
||||
// * With full assembly, use the BoomerAMG preconditioner from hypre.
|
||||
// * With partial assembly, use Jacobi smoothing, for now.
|
||||
Solver *prec = NULL;
|
||||
if (pa_)
|
||||
{
|
||||
if (UsesTensorBasis(fespace))
|
||||
{
|
||||
prec = new OperatorJacobiSmoother(a, ess_tdof_list);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
prec = new HypreBoomerAMG;
|
||||
}
|
||||
CGSolver cg(MPI_COMM_WORLD);
|
||||
cg.SetRelTol(1e-12);
|
||||
cg.SetMaxIter(2000);
|
||||
cg.SetPrintLevel(1);
|
||||
if (prec) { cg.SetPreconditioner(*prec); }
|
||||
cg.SetOperator(*A);
|
||||
cg.Mult(B, X);
|
||||
delete prec;
|
||||
|
||||
// 14. Recover the parallel grid function corresponding to X. This is the
|
||||
// local finite element solution on each processor.
|
||||
a.RecoverFEMSolution(X, b, x);
|
||||
}
|
||||
|
||||
void SetRanges()
|
||||
{
|
||||
switch (coords_)
|
||||
{
|
||||
case POLAR:
|
||||
if (isnan(q1_min_)) { q1_min_ = 0.5; }
|
||||
if (isnan(q1_max_)) { q1_max_ = 4.0; }
|
||||
if (isnan(q2_min_)) { q2_min_ = -M_PI; }
|
||||
if (isnan(q2_max_)) { q2_max_ = M_PI; }
|
||||
break;
|
||||
case PARABOLIC_CYL:
|
||||
if (isnan(q1_min_)) { q1_min_ = -4.0; }
|
||||
if (isnan(q1_max_)) { q1_max_ = 4.0; }
|
||||
if (isnan(q2_min_)) { q2_min_ = 0.5; }
|
||||
if (isnan(q2_max_)) { q2_max_ = 4.0; }
|
||||
break;
|
||||
case ELLIPTIC:
|
||||
if (isnan(q1_min_)) { q1_min_ = 0.5; }
|
||||
if (isnan(q1_max_)) { q1_max_ = 2.0; }
|
||||
if (isnan(q2_min_)) { q2_min_ = -M_PI; }
|
||||
if (isnan(q2_max_)) { q2_max_ = M_PI; }
|
||||
break;
|
||||
case BIPOLAR:
|
||||
if (isnan(q1_min_)) { q1_min_ = -M_PI; }
|
||||
if (isnan(q1_max_)) { q1_max_ = M_PI; }
|
||||
if (isnan(q2_min_)) { q2_min_ = 0.5; }
|
||||
if (isnan(q2_max_)) { q2_max_ = 4.0; }
|
||||
break;
|
||||
case CYLINDRICAL:
|
||||
if (isnan(q1_min_)) { q1_min_ = 0.0; }
|
||||
if (isnan(q1_max_)) { q1_max_ = 4.0; }
|
||||
if (isnan(q2_min_)) { q2_min_ = 0.0; }
|
||||
if (isnan(q2_max_)) { q2_max_ = 4.0; }
|
||||
break;
|
||||
case SPHERICAL:
|
||||
if (isnan(q1_min_)) { q1_min_ = 0.5; }
|
||||
if (isnan(q1_max_)) { q1_max_ = 4.0; }
|
||||
if (isnan(q2_min_)) { q2_min_ = 0.0; }
|
||||
if (isnan(q2_max_)) { q2_max_ = M_PI; }
|
||||
break;
|
||||
case PARABOLIC:
|
||||
if (isnan(q1_min_)) { q1_min_ = 0.2; }
|
||||
if (isnan(q1_max_)) { q1_max_ = 4.0; }
|
||||
if (isnan(q2_min_)) { q2_min_ = 0.2; }
|
||||
if (isnan(q2_max_)) { q2_max_ = 4.0; }
|
||||
break;
|
||||
case PROLATE_SPHEROIDAL:
|
||||
if (isnan(q1_min_)) { q1_min_ = 0.5; }
|
||||
if (isnan(q1_max_)) { q1_max_ = 2.0; }
|
||||
if (isnan(q2_min_)) { q2_min_ = 0.0; }
|
||||
if (isnan(q2_max_)) { q2_max_ = M_PI; }
|
||||
break;
|
||||
case OBLATE_SPHEROIDAL:
|
||||
if (isnan(q1_min_)) { q1_min_ = 0.5; }
|
||||
if (isnan(q1_max_)) { q1_max_ = 2.0; }
|
||||
if (isnan(q2_min_)) { q2_min_ = -0.5 * M_PI; }
|
||||
if (isnan(q2_max_)) { q2_max_ = 0.5 * M_PI; }
|
||||
break;
|
||||
case TOROIDAL:
|
||||
if (isnan(q1_min_)) { q1_min_ = -M_PI; }
|
||||
if (isnan(q1_max_)) { q1_max_ = M_PI; }
|
||||
if (isnan(q2_min_)) { q2_min_ = 0.5; }
|
||||
if (isnan(q2_max_)) { q2_max_ = 4.0; }
|
||||
break;
|
||||
case BISPHERICAL:
|
||||
if (isnan(q1_min_)) { q1_min_ = 0.0; }
|
||||
if (isnan(q1_max_)) { q1_max_ = M_PI; }
|
||||
if (isnan(q2_min_)) { q2_min_ = 0.5; }
|
||||
if (isnan(q2_max_)) { q2_max_ = 4.0; }
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
void trans(const Vector &u, Vector &x)
|
||||
{
|
||||
x.SetSize(2);
|
||||
|
||||
switch (coords_)
|
||||
{
|
||||
case POLAR:
|
||||
x[0] = u[0] * cos(u[1]);
|
||||
x[1] = u[0] * sin(u[1]);
|
||||
break;
|
||||
case PARABOLIC_CYL:
|
||||
x[0] = 0.5 * (u[0] * u[0] - u[1] * u[1]);
|
||||
x[1] = u[0] * u[1];
|
||||
break;
|
||||
case ELLIPTIC:
|
||||
x[0] = a_ * cosh(u[0]) * cos(u[1]);
|
||||
x[1] = a_ * sinh(u[0]) * sin(u[1]);
|
||||
break;
|
||||
case BIPOLAR:
|
||||
{
|
||||
double den = (cosh(u[1]) - cos(u[0]));
|
||||
x[0] = a_ * sinh(u[1]) / den;
|
||||
x[1] = a_ * sin(u[0]) / den;
|
||||
}
|
||||
break;
|
||||
case CYLINDRICAL:
|
||||
{
|
||||
x[0] = u[0];
|
||||
x[1] = u[1];
|
||||
}
|
||||
break;
|
||||
case SPHERICAL:
|
||||
{
|
||||
x[0] = u[0] * sin(u[1]);
|
||||
x[1] = -u[0] * cos(u[1]);
|
||||
}
|
||||
break;
|
||||
case PARABOLIC:
|
||||
{
|
||||
x[0] = u[0] * u[1];
|
||||
x[1] = -0.5 * (u[0] * u[0] - u[1] * u[1]);
|
||||
}
|
||||
break;
|
||||
case PROLATE_SPHEROIDAL:
|
||||
{
|
||||
x[0] = a_ * sinh(u[0]) * sin(u[1]);
|
||||
x[1] = -a_ * cosh(u[0]) * cos(u[1]);
|
||||
}
|
||||
break;
|
||||
case OBLATE_SPHEROIDAL:
|
||||
{
|
||||
x[0] = a_ * cosh(u[0]) * cos(u[1]);
|
||||
x[1] = a_ * sinh(u[0]) * sin(u[1]);
|
||||
}
|
||||
break;
|
||||
case TOROIDAL:
|
||||
{
|
||||
double den = (cosh(u[1]) - cos(u[0]));
|
||||
x[0] = a_ * sinh(u[1]) / den;
|
||||
x[1] = a_ * sin(u[0]) / den;
|
||||
}
|
||||
break;
|
||||
case BISPHERICAL:
|
||||
{
|
||||
double den = (cosh(u[1]) + cos(u[0]));
|
||||
x[0] = a_ * sin(u[0]) / den;
|
||||
x[1] = a_ * sinh(u[1]) / den;
|
||||
}
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
double OrthoCoef::Eval(ElementTransformation &T,
|
||||
const IntegrationPoint &ip)
|
||||
{
|
||||
double u_data[2];
|
||||
Vector u(u_data, 2);
|
||||
T.Transform(ip, u);
|
||||
|
||||
switch (coords_)
|
||||
{
|
||||
case POLAR:
|
||||
switch (ind_)
|
||||
{
|
||||
case 0:
|
||||
return 1.0;
|
||||
case 1:
|
||||
return u[0];
|
||||
case 2:
|
||||
return 1.0;
|
||||
default:
|
||||
return 0.0;
|
||||
}
|
||||
break;
|
||||
case PARABOLIC_CYL:
|
||||
switch (ind_)
|
||||
{
|
||||
case 0:
|
||||
return sqrt(u[0] * u[0] + u[1] * u[1]);
|
||||
case 1:
|
||||
return sqrt(u[0] * u[0] + u[1] * u[1]);
|
||||
case 2:
|
||||
return 1.0;
|
||||
default:
|
||||
return 0.0;
|
||||
}
|
||||
break;
|
||||
case ELLIPTIC:
|
||||
switch (ind_)
|
||||
{
|
||||
case 0:
|
||||
return a_ * sqrt(pow(sinh(u[0]), 2) + pow(sin(u[1]), 2));
|
||||
case 1:
|
||||
return a_ * sqrt(pow(sinh(u[0]), 2) + pow(sin(u[1]), 2));
|
||||
case 2:
|
||||
return 1.0;
|
||||
default:
|
||||
return 0.0;
|
||||
}
|
||||
break;
|
||||
case BIPOLAR:
|
||||
{
|
||||
double den = (cosh(u[1]) - cos(u[0]));
|
||||
|
||||
switch (ind_)
|
||||
{
|
||||
case 0:
|
||||
return a_ / den;
|
||||
case 1:
|
||||
return a_ / den;
|
||||
case 2:
|
||||
return 1.0;
|
||||
default:
|
||||
return 0.0;
|
||||
}
|
||||
break;
|
||||
}
|
||||
case CYLINDRICAL:
|
||||
switch (ind_)
|
||||
{
|
||||
case 0:
|
||||
return 1.0;
|
||||
case 1:
|
||||
return 1.0;
|
||||
case 2:
|
||||
return u[0];
|
||||
default:
|
||||
return 0.0;
|
||||
}
|
||||
break;
|
||||
case SPHERICAL:
|
||||
switch (ind_)
|
||||
{
|
||||
case 0:
|
||||
return 1.0;
|
||||
case 1:
|
||||
return u[0];
|
||||
case 2:
|
||||
return u[0] * sin(u[1]);
|
||||
default:
|
||||
return 0.0;
|
||||
}
|
||||
break;
|
||||
case PARABOLIC:
|
||||
switch (ind_)
|
||||
{
|
||||
case 0:
|
||||
return sqrt(u[0] * u[0] + u[1] * u[1]);
|
||||
case 1:
|
||||
return sqrt(u[0] * u[0] + u[1] * u[1]);
|
||||
case 2:
|
||||
return u[0] * u[1];
|
||||
default:
|
||||
return 0.0;
|
||||
}
|
||||
break;
|
||||
case PROLATE_SPHEROIDAL:
|
||||
switch (ind_)
|
||||
{
|
||||
case 0:
|
||||
return a_ * sqrt(pow(sinh(u[0]), 2) + pow(sin(u[1]), 2));
|
||||
case 1:
|
||||
return a_ * sqrt(pow(sinh(u[0]), 2) + pow(sin(u[1]), 2));
|
||||
case 2:
|
||||
return a_ * sinh(u[0]) * sin(u[1]);
|
||||
default:
|
||||
return 0.0;
|
||||
}
|
||||
break;
|
||||
case OBLATE_SPHEROIDAL:
|
||||
switch (ind_)
|
||||
{
|
||||
case 0:
|
||||
return a_ * sqrt(pow(sinh(u[0]), 2) + pow(sin(u[1]), 2));
|
||||
case 1:
|
||||
return a_ * sqrt(pow(sinh(u[0]), 2) + pow(sin(u[1]), 2));
|
||||
case 2:
|
||||
return a_ * cosh(u[0]) * cos(u[1]);
|
||||
default:
|
||||
return 0.0;
|
||||
}
|
||||
break;
|
||||
case TOROIDAL:
|
||||
{
|
||||
double den = (cosh(u[1]) - cos(u[0]));
|
||||
|
||||
switch (ind_)
|
||||
{
|
||||
case 0:
|
||||
return a_ / den;
|
||||
case 1:
|
||||
return a_ / den;
|
||||
case 2:
|
||||
return a_ * sinh(u[1]) / den;
|
||||
default:
|
||||
return 0.0;
|
||||
}
|
||||
break;
|
||||
}
|
||||
case BISPHERICAL:
|
||||
{
|
||||
double den = (cosh(u[1]) + cos(u[0]));
|
||||
|
||||
switch (ind_)
|
||||
{
|
||||
case 0:
|
||||
return a_ / den;
|
||||
case 1:
|
||||
return a_ / den;
|
||||
case 2:
|
||||
return a_ * sin(u[0]) / den;
|
||||
default:
|
||||
return 0.0;
|
||||
}
|
||||
break;
|
||||
}
|
||||
}
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
double OrthoWeightCoef::Eval(ElementTransformation &T,
|
||||
const IntegrationPoint &ip)
|
||||
{
|
||||
double h1 = h1Coef_.Eval(T, ip);
|
||||
double h2 = h2Coef_.Eval(T, ip);
|
||||
double h3 = h3Coef_.Eval(T, ip);
|
||||
|
||||
return h1 * h2 * h3;
|
||||
}
|
||||
|
||||
void OrthoMatrixCoef::Eval(DenseMatrix &K, ElementTransformation &T,
|
||||
const IntegrationPoint &ip)
|
||||
{
|
||||
double h1 = h1Coef_.Eval(T, ip);
|
||||
double h2 = h2Coef_.Eval(T, ip);
|
||||
double h3 = h3Coef_.Eval(T, ip);
|
||||
|
||||
K.SetSize(2);
|
||||
K(0,0) = h2 * h3 / h1;
|
||||
K(0,1) = 0.0;
|
||||
K(1,0) = 0.0;
|
||||
K(1,1) = h1 * h3 / h2;
|
||||
}
|
||||
+1
-1
@@ -26,7 +26,7 @@ SEQ_EXAMPLES = ex0 ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 \
|
||||
PAR_EXAMPLES = ex0p ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p \
|
||||
ex12p ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p \
|
||||
ex25p ex26p ex27p ex28p ex29p ex30p ex31p ex32p ex33p ex34p ex35p ex36p \
|
||||
ex37p ex39p ex40p ex1p-orth ex11p-cyl ex13p-cyl ex13p-cyl-3d
|
||||
ex37p ex39p ex40p
|
||||
SEQ_DEVICE_EXAMPLES = ex1 ex3 ex4 ex5 ex6 ex9 ex14 ex22 ex24 ex25 ex26 ex34
|
||||
PAR_DEVICE_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex14p \
|
||||
ex22p ex24p ex25p ex26p ex34p ex35p
|
||||
|
||||
@@ -1,214 +0,0 @@
|
||||
// Copyright (c) 2010-2025, 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.
|
||||
|
||||
|
||||
// Abstract array data type
|
||||
|
||||
#include "array.hpp"
|
||||
#include "../general/forall.hpp"
|
||||
#include <fstream>
|
||||
#include <type_traits>
|
||||
|
||||
namespace mfem
|
||||
{
|
||||
|
||||
template <class T>
|
||||
void Array<T>::Print(std::ostream &os, int width) const
|
||||
{
|
||||
for (int i = 0; i < size; i++)
|
||||
{
|
||||
os << data[i];
|
||||
if ( !((i+1) % width) || i+1 == size )
|
||||
{
|
||||
os << '\n';
|
||||
}
|
||||
else
|
||||
{
|
||||
os << " ";
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <class T>
|
||||
void Array<T>::Save(std::ostream &os, int fmt) const
|
||||
{
|
||||
if (fmt == 0)
|
||||
{
|
||||
os << size << '\n';
|
||||
}
|
||||
for (int i = 0; i < size; i++)
|
||||
{
|
||||
os << operator[](i) << '\n';
|
||||
}
|
||||
}
|
||||
|
||||
template <class T>
|
||||
void Array<T>::Load(std::istream &in, int fmt)
|
||||
{
|
||||
if (fmt == 0)
|
||||
{
|
||||
int new_size;
|
||||
in >> new_size;
|
||||
SetSize(new_size);
|
||||
}
|
||||
for (int i = 0; i < size; i++)
|
||||
{
|
||||
in >> operator[](i);
|
||||
}
|
||||
}
|
||||
|
||||
template <class T>
|
||||
T Array<T>::Max() const
|
||||
{
|
||||
MFEM_ASSERT(size > 0, "Array is empty with size " << size);
|
||||
|
||||
T max = operator[](0);
|
||||
for (int i = 1; i < size; i++)
|
||||
{
|
||||
if (max < operator[](i))
|
||||
{
|
||||
max = operator[](i);
|
||||
}
|
||||
}
|
||||
|
||||
return max;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
T Array<T>::Min() const
|
||||
{
|
||||
MFEM_ASSERT(size > 0, "Array is empty with size " << size);
|
||||
|
||||
T min = operator[](0);
|
||||
for (int i = 1; i < size; i++)
|
||||
{
|
||||
if (operator[](i) < min)
|
||||
{
|
||||
min = operator[](i);
|
||||
}
|
||||
}
|
||||
|
||||
return min;
|
||||
}
|
||||
|
||||
// Partial Sum
|
||||
template <class T>
|
||||
void Array<T>::PartialSum()
|
||||
{
|
||||
T sum = static_cast<T>(0);
|
||||
for (int i = 0; i < size; i++)
|
||||
{
|
||||
sum+=operator[](i);
|
||||
operator[](i) = sum;
|
||||
}
|
||||
}
|
||||
|
||||
template <class T>
|
||||
void Array<T>::Abs()
|
||||
{
|
||||
static_assert(std::is_arithmetic<T>::value, "Use with arithmetic types!");
|
||||
const bool useDevice = UseDevice();
|
||||
const int N = size;
|
||||
auto y = ReadWrite(useDevice);
|
||||
mfem::forall_switch(useDevice, N, [=] MFEM_HOST_DEVICE (int i)
|
||||
{
|
||||
y[i] = std::abs(y[i]);
|
||||
});
|
||||
}
|
||||
|
||||
// Sum
|
||||
template <class T>
|
||||
T Array<T>::Sum() const
|
||||
{
|
||||
T sum = static_cast<T>(0);
|
||||
for (int i = 0; i < size; i++)
|
||||
{
|
||||
sum+=operator[](i);
|
||||
}
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
int Array<T>::IsSorted() const
|
||||
{
|
||||
T val_prev = operator[](0), val;
|
||||
for (int i = 1; i < size; i++)
|
||||
{
|
||||
val=operator[](i);
|
||||
if (val < val_prev)
|
||||
{
|
||||
return 0;
|
||||
}
|
||||
val_prev = val;
|
||||
}
|
||||
|
||||
return 1;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
bool Array<T>::IsConstant() const
|
||||
{
|
||||
if (size < 2) { return true; }
|
||||
const T v0 = data[0];
|
||||
for (int i = 1; i < size; i++)
|
||||
{
|
||||
if (data[i] != v0)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
void Array2D<T>::Load(const char *filename, int fmt)
|
||||
{
|
||||
std::ifstream in;
|
||||
in.open(filename, std::ifstream::in);
|
||||
MFEM_VERIFY(in.is_open(), "File " << filename << " does not exist.");
|
||||
Load(in, fmt);
|
||||
in.close();
|
||||
}
|
||||
|
||||
template <class T>
|
||||
void Array2D<T>::Print(std::ostream &os, int width_)
|
||||
{
|
||||
int height = this->NumRows();
|
||||
int width = this->NumCols();
|
||||
|
||||
for (int i = 0; i < height; i++)
|
||||
{
|
||||
os << "[row " << i << "]\n";
|
||||
for (int j = 0; j < width; j++)
|
||||
{
|
||||
os << (*this)(i,j);
|
||||
if ( (j+1) == width_ || (j+1) % width_ == 0 )
|
||||
{
|
||||
os << '\n';
|
||||
}
|
||||
else
|
||||
{
|
||||
os << ' ';
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template class Array<char>;
|
||||
template class Array<int>;
|
||||
template class Array<long long>;
|
||||
template class Array<real_t>;
|
||||
template class Array2D<int>;
|
||||
template class Array2D<real_t>;
|
||||
|
||||
} // namespace mfem
|
||||
+213
-15
@@ -16,9 +16,13 @@
|
||||
#include "mem_manager.hpp"
|
||||
#include "device.hpp"
|
||||
#include "error.hpp"
|
||||
#include "forall.hpp"
|
||||
#include "globals.hpp"
|
||||
#include "reducers.hpp"
|
||||
#include "scan.hpp"
|
||||
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <cstdlib>
|
||||
#include <cstring>
|
||||
#include <algorithm>
|
||||
@@ -135,6 +139,8 @@ public:
|
||||
/// Return the device flag of the Memory object used by the Array
|
||||
bool UseDevice() const { return data.UseDevice(); }
|
||||
|
||||
void UseDevice(bool use_dev) { data.UseDevice(use_dev); }
|
||||
|
||||
/// Return true if the data will be deleted by the Array
|
||||
inline bool OwnsData() const { return data.OwnsHostPtr(); }
|
||||
|
||||
@@ -275,11 +281,11 @@ public:
|
||||
|
||||
/** @brief Find the maximal element in the array, using the comparison
|
||||
operator `<` for class T. */
|
||||
T Max() const;
|
||||
inline T Max() const;
|
||||
|
||||
/** @brief Find the minimal element in the array, using the comparison
|
||||
operator `<` for class T. */
|
||||
T Min() const;
|
||||
inline T Min() const;
|
||||
|
||||
/// Sorts the array in ascending order. This requires operator< to be defined for T.
|
||||
void Sort() { std::sort((T*)data, data + size); }
|
||||
@@ -297,22 +303,22 @@ public:
|
||||
}
|
||||
|
||||
/// Return 1 if the array is sorted from lowest to highest. Otherwise return 0.
|
||||
int IsSorted() const;
|
||||
inline int IsSorted() const;
|
||||
|
||||
/// Does the Array have Size zero.
|
||||
bool IsEmpty() const { return Size() == 0; }
|
||||
|
||||
/// Return true if all entries of the array are the same.
|
||||
bool IsConstant() const;
|
||||
inline bool IsConstant() const;
|
||||
|
||||
/// Fill the entries of the array with the cumulative sum of the entries.
|
||||
void PartialSum();
|
||||
inline void PartialSum();
|
||||
|
||||
/// Replace each entry of the array with its absolute value.
|
||||
void Abs();
|
||||
inline void Abs();
|
||||
|
||||
/// Return the sum of all the array entries using the '+'' operator for class 'T'.
|
||||
T Sum() const;
|
||||
inline T Sum() const;
|
||||
|
||||
/// Set all entries of the array to the provided constant.
|
||||
inline void operator=(const T &a);
|
||||
@@ -797,8 +803,14 @@ template <typename T> template <typename CT>
|
||||
inline Array<T> &Array<T>::operator=(const Array<CT> &src)
|
||||
{
|
||||
SetSize(src.Size());
|
||||
for (int i = 0; i < size; i++) { (*this)[i] = T(src[i]); }
|
||||
return *this;
|
||||
|
||||
const bool use_dev = UseDevice() || src.UseDevice();
|
||||
const auto x = src.Read(use_dev);
|
||||
auto y = Write(use_dev);
|
||||
mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
|
||||
{
|
||||
y[i] = x[i];
|
||||
});
|
||||
}
|
||||
|
||||
template <class T>
|
||||
@@ -1014,19 +1026,24 @@ template <class T>
|
||||
inline void Array<T>::GetSubArray(int offset, int sa_size, Array<T> &sa) const
|
||||
{
|
||||
sa.SetSize(sa_size);
|
||||
for (int i = 0; i < sa_size; i++)
|
||||
const bool use_dev = UseDevice() || sa.UseDevice();
|
||||
const auto x = Read(use_dev);
|
||||
auto y = sa.Write(use_dev);
|
||||
mfem::forall_switch(use_dev, sa_size, [=] MFEM_HOST_DEVICE (int i)
|
||||
{
|
||||
sa[i] = (*this)[offset+i];
|
||||
}
|
||||
y[i] = x[offset + i];
|
||||
});
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline void Array<T>::operator=(const T &a)
|
||||
{
|
||||
for (int i = 0; i < size; i++)
|
||||
const bool use_dev = UseDevice();
|
||||
auto x = Write(use_dev);
|
||||
mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
|
||||
{
|
||||
data[i] = a;
|
||||
}
|
||||
x[i] = a;
|
||||
});
|
||||
}
|
||||
|
||||
template <class T>
|
||||
@@ -1035,6 +1052,153 @@ inline void Array<T>::Assign(const T *p)
|
||||
data.CopyFromHost(p, Size());
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline void Array<T>::Print(std::ostream &os, int width) const
|
||||
{
|
||||
for (int i = 0; i < size; i++)
|
||||
{
|
||||
os << data[i];
|
||||
if ( !((i+1) % width) || i+1 == size )
|
||||
{
|
||||
os << '\n';
|
||||
}
|
||||
else
|
||||
{
|
||||
os << " ";
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline void Array<T>::Save(std::ostream &os, int fmt) const
|
||||
{
|
||||
if (fmt == 0)
|
||||
{
|
||||
os << size << '\n';
|
||||
}
|
||||
for (int i = 0; i < size; i++)
|
||||
{
|
||||
os << operator[](i) << '\n';
|
||||
}
|
||||
}
|
||||
|
||||
template <class T>
|
||||
void Array<T>::Load(std::istream &in, int fmt)
|
||||
{
|
||||
if (fmt == 0)
|
||||
{
|
||||
int new_size;
|
||||
in >> new_size;
|
||||
SetSize(new_size);
|
||||
}
|
||||
for (int i = 0; i < size; i++)
|
||||
{
|
||||
in >> operator[](i);
|
||||
}
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline T Array<T>::Max() const
|
||||
{
|
||||
MFEM_ASSERT(size > 0, "Array is empty with size " << size);
|
||||
|
||||
T max = operator[](0);
|
||||
for (int i = 1; i < size; i++)
|
||||
{
|
||||
if (max < operator[](i))
|
||||
{
|
||||
max = operator[](i);
|
||||
}
|
||||
}
|
||||
|
||||
return max;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline T Array<T>::Min() const
|
||||
{
|
||||
MFEM_ASSERT(size > 0, "Array is empty with size " << size);
|
||||
|
||||
T min = operator[](0);
|
||||
for (int i = 1; i < size; i++)
|
||||
{
|
||||
if (operator[](i) < min)
|
||||
{
|
||||
min = operator[](i);
|
||||
}
|
||||
}
|
||||
|
||||
return min;
|
||||
}
|
||||
|
||||
// Partial Sum
|
||||
template <class T>
|
||||
inline void Array<T>::PartialSum()
|
||||
{
|
||||
auto data_ptr = ReadWrite(UseDevice());
|
||||
InclusiveScan(UseDevice(), data_ptr, data_ptr, size);
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline void Array<T>::Abs()
|
||||
{
|
||||
static_assert(std::is_arithmetic<T>::value, "Use with arithmetic types!");
|
||||
const bool useDevice = UseDevice();
|
||||
const int N = size;
|
||||
auto y = ReadWrite(useDevice);
|
||||
mfem::forall_switch(useDevice, N, [=] MFEM_HOST_DEVICE (int i)
|
||||
{
|
||||
y[i] = std::abs(y[i]);
|
||||
});
|
||||
}
|
||||
|
||||
// Sum
|
||||
template <class T>
|
||||
inline T Array<T>::Sum() const
|
||||
{
|
||||
T sum = static_cast<T>(0);
|
||||
if (size > 0)
|
||||
{
|
||||
const auto m_data = Read(UseDevice());
|
||||
reduce(size, sum, [=] MFEM_HOST_DEVICE(int i, T &r) { r += m_data[i]; },
|
||||
/* */ SumReducer<T> {}, UseDevice());
|
||||
}
|
||||
return sum;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline int Array<T>::IsSorted() const
|
||||
{
|
||||
T val_prev = operator[](0), val;
|
||||
for (int i = 1; i < size; i++)
|
||||
{
|
||||
val=operator[](i);
|
||||
if (val < val_prev)
|
||||
{
|
||||
return 0;
|
||||
}
|
||||
val_prev = val;
|
||||
}
|
||||
|
||||
return 1;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline bool Array<T>::IsConstant() const
|
||||
{
|
||||
if (size < 2) { return true; }
|
||||
const T v0 = data[0];
|
||||
for (int i = 1; i < size; i++)
|
||||
{
|
||||
if (data[i] != v0)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
template <class T>
|
||||
inline const T &Array2D<T>::operator()(int i, int j) const
|
||||
@@ -1074,6 +1238,40 @@ inline T *Array2D<T>::operator[](int i)
|
||||
return &array1d[i*N];
|
||||
}
|
||||
|
||||
template <class T>
|
||||
void Array2D<T>::Load(const char *filename, int fmt)
|
||||
{
|
||||
std::ifstream in;
|
||||
in.open(filename, std::ifstream::in);
|
||||
MFEM_VERIFY(in.is_open(), "File " << filename << " does not exist.");
|
||||
Load(in, fmt);
|
||||
in.close();
|
||||
}
|
||||
|
||||
template <class T>
|
||||
void Array2D<T>::Print(std::ostream &os, int width_)
|
||||
{
|
||||
int height = this->NumRows();
|
||||
int width = this->NumCols();
|
||||
|
||||
for (int i = 0; i < height; i++)
|
||||
{
|
||||
os << "[row " << i << "]\n";
|
||||
for (int j = 0; j < width; j++)
|
||||
{
|
||||
os << (*this)(i,j);
|
||||
if ( (j+1) == width_ || (j+1) % width_ == 0 )
|
||||
{
|
||||
os << '\n';
|
||||
}
|
||||
else
|
||||
{
|
||||
os << ' ';
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template <class T>
|
||||
inline void Swap(Array2D<T> &a, Array2D<T> &b)
|
||||
|
||||
+29
-10
@@ -12,7 +12,6 @@
|
||||
#ifndef MFEM_REDUCERS_HPP
|
||||
#define MFEM_REDUCERS_HPP
|
||||
|
||||
#include "array.hpp"
|
||||
#include "forall.hpp"
|
||||
|
||||
#include <cmath>
|
||||
@@ -514,6 +513,33 @@ template<class B, class R> struct reduction_kernel
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
template <class T>
|
||||
class ReductionWorkspace
|
||||
{
|
||||
Memory<T> workspace;
|
||||
|
||||
static ReductionWorkspace &Instance()
|
||||
{
|
||||
static ReductionWorkspace instance;
|
||||
return instance;
|
||||
}
|
||||
|
||||
~ReductionWorkspace() { workspace.Delete(); }
|
||||
|
||||
public:
|
||||
static T *Get(int num_blocks)
|
||||
{
|
||||
ReductionWorkspace &instance = Instance();
|
||||
if (instance.workspace.Capacity() < num_blocks)
|
||||
{
|
||||
instance.workspace.Delete();
|
||||
instance.workspace.New(num_blocks, MemoryType::HOST_PINNED);
|
||||
}
|
||||
return instance.workspace;
|
||||
}
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
@@ -529,8 +555,7 @@ template<class B, class R> struct reduction_kernel
|
||||
@tparam T value_type to operate on
|
||||
*/
|
||||
template <class T, class B, class R>
|
||||
void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,
|
||||
Array<T> &workspace)
|
||||
void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev)
|
||||
{
|
||||
if (N == 0)
|
||||
{
|
||||
@@ -567,13 +592,7 @@ void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,
|
||||
|
||||
red_type red{nullptr, std::forward<B>(body), reducer, N, items_per_thread};
|
||||
// allocate res to fit block_size entries
|
||||
auto mt = workspace.GetMemory().GetMemoryType();
|
||||
if (mt != MemoryType::HOST_PINNED && mt != MemoryType::MANAGED)
|
||||
{
|
||||
mt = MemoryType::HOST_PINNED;
|
||||
}
|
||||
workspace.SetSize(nblocks, mt);
|
||||
auto work = workspace.HostWrite();
|
||||
auto work = internal::ReductionWorkspace<T>::Get(nblocks);
|
||||
red.work = work;
|
||||
forall_2D(nblocks, block_size, 1, std::move(red));
|
||||
// wait for results
|
||||
|
||||
+52
-22
@@ -28,8 +28,37 @@
|
||||
|
||||
namespace mfem
|
||||
{
|
||||
/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, workspace,
|
||||
/// std::plus<>{})
|
||||
|
||||
namespace internal
|
||||
{
|
||||
class ScanWorkspace
|
||||
{
|
||||
Memory<std::byte> workspace;
|
||||
static ScanWorkspace &Instance()
|
||||
{
|
||||
static ScanWorkspace instance;
|
||||
return instance;
|
||||
}
|
||||
~ScanWorkspace() { workspace.Delete(); }
|
||||
public:
|
||||
static std::byte *Get(int num_bytes)
|
||||
{
|
||||
ScanWorkspace &instance = Instance();
|
||||
if (Size() < num_bytes)
|
||||
{
|
||||
instance.workspace.Delete();
|
||||
instance.workspace.New(num_bytes);
|
||||
}
|
||||
return instance.workspace.Write(MemoryClass::DEVICE, Size());
|
||||
}
|
||||
static int Size()
|
||||
{
|
||||
return Instance().workspace.Capacity();
|
||||
}
|
||||
};
|
||||
}
|
||||
|
||||
/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, std::plus<>{})
|
||||
template <class InputIt, class OutputIt>
|
||||
void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
|
||||
{
|
||||
@@ -37,12 +66,12 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
|
||||
#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
|
||||
if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
|
||||
{
|
||||
static Array<std::byte> workspace;
|
||||
size_t bytes = workspace.Size();
|
||||
if (bytes)
|
||||
using internal::ScanWorkspace;
|
||||
size_t bytes = ScanWorkspace::Size();
|
||||
if (bytes > 0)
|
||||
{
|
||||
auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
|
||||
workspace.Write(), bytes, d_in, d_out, num_items);
|
||||
ScanWorkspace::Get(bytes), bytes, d_in, d_out, num_items);
|
||||
#if defined(MFEM_USE_CUDA)
|
||||
if (err == cudaSuccess)
|
||||
{
|
||||
@@ -57,11 +86,12 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
|
||||
}
|
||||
// try allocating a larger buffer
|
||||
bytes = 0;
|
||||
// get size of buffer
|
||||
MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
|
||||
nullptr, bytes, d_in, d_out, num_items));
|
||||
workspace.SetSize(bytes);
|
||||
// resize buffer (in ScanWorkspace::Get) and try again
|
||||
MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
|
||||
workspace.Write(), bytes, d_in, d_out, num_items));
|
||||
ScanWorkspace::Get(bytes), bytes, d_in, d_out, num_items));
|
||||
return;
|
||||
}
|
||||
#endif
|
||||
@@ -101,12 +131,13 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
|
||||
#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
|
||||
if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
|
||||
{
|
||||
static Array<std::byte> workspace;
|
||||
size_t bytes = workspace.Size();
|
||||
if (bytes)
|
||||
using internal::ScanWorkspace;
|
||||
size_t bytes = ScanWorkspace::Size();
|
||||
if (bytes > 0)
|
||||
{
|
||||
auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
|
||||
workspace.Write(), bytes, d_in, d_out, scan_op, num_items);
|
||||
ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
|
||||
num_items);
|
||||
#if defined(MFEM_USE_CUDA)
|
||||
if (err == cudaSuccess)
|
||||
{
|
||||
@@ -123,9 +154,9 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
|
||||
bytes = 0;
|
||||
MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
|
||||
nullptr, bytes, d_in, d_out, scan_op, num_items));
|
||||
workspace.SetSize(bytes);
|
||||
MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
|
||||
workspace.Write(), bytes, d_in, d_out, scan_op, num_items));
|
||||
ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
|
||||
num_items));
|
||||
return;
|
||||
}
|
||||
#endif
|
||||
@@ -164,13 +195,13 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
|
||||
#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
|
||||
if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
|
||||
{
|
||||
static Array<std::byte> workspace;
|
||||
size_t bytes = workspace.Size();
|
||||
using internal::ScanWorkspace;
|
||||
size_t bytes = ScanWorkspace::Size();
|
||||
if (bytes)
|
||||
{
|
||||
auto err = MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
|
||||
workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
|
||||
num_items);
|
||||
ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
|
||||
init_value, num_items);
|
||||
#if defined(MFEM_USE_CUDA)
|
||||
if (err == cudaSuccess)
|
||||
{
|
||||
@@ -187,10 +218,9 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
|
||||
bytes = 0;
|
||||
MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
|
||||
nullptr, bytes, d_in, d_out, scan_op, init_value, num_items));
|
||||
workspace.SetSize(bytes);
|
||||
MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
|
||||
workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
|
||||
num_items));
|
||||
ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
|
||||
init_value, num_items));
|
||||
return;
|
||||
}
|
||||
#endif
|
||||
@@ -213,7 +243,7 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
|
||||
}
|
||||
|
||||
/// Equivalent to ExclusiveScan(use_dev, d_in, d_out, num_items, init_value,
|
||||
/// workspace, std::plus<>{})
|
||||
/// std::plus<>{})
|
||||
template <class InputIt, class OutputIt, class T>
|
||||
void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
|
||||
T init_value)
|
||||
|
||||
+8
-20
@@ -92,18 +92,6 @@ struct LpReducer
|
||||
}
|
||||
};
|
||||
|
||||
static Array<real_t>& vector_workspace()
|
||||
{
|
||||
static Array<real_t> instance;
|
||||
return instance;
|
||||
}
|
||||
|
||||
static Array<DevicePair<real_t, real_t>> &Lpvector_workspace()
|
||||
{
|
||||
static Array<DevicePair<real_t, real_t>> instance;
|
||||
return instance;
|
||||
}
|
||||
|
||||
Vector::Vector(const Vector &v)
|
||||
{
|
||||
const int s = v.Size();
|
||||
@@ -991,7 +979,7 @@ real_t Vector::Norml2() const
|
||||
}
|
||||
}
|
||||
},
|
||||
L2Reducer{}, UseDevice(), Lpvector_workspace());
|
||||
L2Reducer{}, UseDevice());
|
||||
// final answer
|
||||
return res.second * sqrt(res.first);
|
||||
}
|
||||
@@ -1006,7 +994,7 @@ real_t Vector::Normlinf() const
|
||||
{
|
||||
r = fmax(r, fabs(m_data[i]));
|
||||
},
|
||||
MaxReducer<real_t> {}, UseDevice(), vector_workspace());
|
||||
MaxReducer<real_t> {}, UseDevice());
|
||||
return res;
|
||||
}
|
||||
|
||||
@@ -1020,7 +1008,7 @@ real_t Vector::Norml1() const
|
||||
{
|
||||
r += fabs(m_data[i]);
|
||||
},
|
||||
SumReducer<real_t> {}, UseDevice(), vector_workspace());
|
||||
SumReducer<real_t> {}, UseDevice());
|
||||
return res;
|
||||
}
|
||||
|
||||
@@ -1063,7 +1051,7 @@ real_t Vector::Normlp(real_t p) const
|
||||
}
|
||||
}
|
||||
},
|
||||
LpReducer{p}, UseDevice(), Lpvector_workspace());
|
||||
LpReducer{p}, UseDevice());
|
||||
// final answer
|
||||
return res.second * pow(res.first, 1.0 / p);
|
||||
} // end if p < infinity()
|
||||
@@ -1096,7 +1084,7 @@ real_t Vector::operator*(const Vector &v) const
|
||||
{
|
||||
r += m_data[i] * v_data[i];
|
||||
},
|
||||
SumReducer<real_t> {}, use_dev, vector_workspace());
|
||||
SumReducer<real_t> {}, use_dev);
|
||||
return res;
|
||||
};
|
||||
|
||||
@@ -1167,7 +1155,7 @@ real_t Vector::Min() const
|
||||
{
|
||||
r = fmin(r, m_data[i]);
|
||||
},
|
||||
MinReducer<real_t> {}, use_dev, vector_workspace());
|
||||
MinReducer<real_t> {}, use_dev);
|
||||
return res;
|
||||
};
|
||||
|
||||
@@ -1213,7 +1201,7 @@ real_t Vector::Max() const
|
||||
{
|
||||
r = fmax(r, m_data[i]);
|
||||
},
|
||||
MaxReducer<real_t> {}, use_dev, vector_workspace());
|
||||
MaxReducer<real_t> {}, use_dev);
|
||||
return res;
|
||||
};
|
||||
|
||||
@@ -1248,7 +1236,7 @@ real_t Vector::Sum() const
|
||||
{
|
||||
r += m_data[i];
|
||||
},
|
||||
SumReducer<real_t> {}, UseDevice(), vector_workspace());
|
||||
SumReducer<real_t> {}, UseDevice());
|
||||
return res;
|
||||
}
|
||||
|
||||
|
||||
@@ -22,7 +22,6 @@ using namespace mfem;
|
||||
|
||||
TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
|
||||
{
|
||||
Array<int> workspace;
|
||||
Array<int> a(1000);
|
||||
a.HostReadWrite();
|
||||
for (int i = 0; i < a.Size(); ++i)
|
||||
@@ -36,7 +35,7 @@ TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
|
||||
int res = 0;
|
||||
mfem::reduce(
|
||||
a.Size(), res, [=] MFEM_HOST_DEVICE(int i, int &r) { r += dptr[i]; },
|
||||
SumReducer<int> {}, use_dev, workspace);
|
||||
SumReducer<int> {}, use_dev);
|
||||
// correct for even-length summations
|
||||
int expected = (AsConst(a)[0] + AsConst(a)[a.Size() - 1]) * a.Size() / 2;
|
||||
CAPTURE(use_dev);
|
||||
@@ -46,7 +45,6 @@ TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
|
||||
|
||||
TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
|
||||
{
|
||||
Array<long long> workspace;
|
||||
Array<long long> a(64);
|
||||
a.HostReadWrite();
|
||||
for (int i = 0; i < a.Size(); ++i)
|
||||
@@ -64,7 +62,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
|
||||
mfem::reduce(
|
||||
a.Size(), res,
|
||||
[=] MFEM_HOST_DEVICE(int i, long long &r) { r *= dptr[i]; },
|
||||
MultReducer<long long> {}, use_dev, workspace);
|
||||
MultReducer<long long> {}, use_dev);
|
||||
long long expected = 0;
|
||||
CAPTURE(use_dev);
|
||||
REQUIRE(res == expected);
|
||||
@@ -76,7 +74,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
|
||||
mfem::reduce(
|
||||
a.Size(), res,
|
||||
[=] MFEM_HOST_DEVICE(int i, long long &r) { r *= dptr[i]; },
|
||||
MultReducer<long long> {}, use_dev, workspace);
|
||||
MultReducer<long long> {}, use_dev);
|
||||
long long expected = 21936950640377856;
|
||||
CAPTURE(use_dev);
|
||||
REQUIRE(res == expected);
|
||||
@@ -86,7 +84,6 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
|
||||
|
||||
TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
|
||||
{
|
||||
Array<unsigned> workspace;
|
||||
Array<unsigned> a(10);
|
||||
SECTION("{ Bit unset }")
|
||||
{
|
||||
@@ -108,7 +105,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
|
||||
mfem::reduce(
|
||||
a.Size(), res,
|
||||
[=] MFEM_HOST_DEVICE(int i, unsigned &r) { r &= dptr[i]; },
|
||||
BAndReducer<unsigned> {}, use_dev, workspace);
|
||||
BAndReducer<unsigned> {}, use_dev);
|
||||
CAPTURE(use_dev);
|
||||
REQUIRE(res == ((~1u) & ~(1u << unset_bit)));
|
||||
REQUIRE((res & (1u << unset_bit)) == 0);
|
||||
@@ -132,7 +129,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
|
||||
mfem::reduce(
|
||||
a.Size(), res,
|
||||
[=] MFEM_HOST_DEVICE(int i, unsigned &r) { r &= dptr[i]; },
|
||||
BAndReducer<unsigned> {}, use_dev, workspace);
|
||||
BAndReducer<unsigned> {}, use_dev);
|
||||
CAPTURE(use_dev);
|
||||
REQUIRE(res == (1u << set_bit));
|
||||
}
|
||||
@@ -141,7 +138,6 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
|
||||
|
||||
TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
|
||||
{
|
||||
Array<unsigned> workspace;
|
||||
Array<unsigned> a(0x210);
|
||||
a.HostReadWrite();
|
||||
for (int i = 0; i < a.Size(); ++i)
|
||||
@@ -157,7 +153,7 @@ TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
|
||||
mfem::reduce(
|
||||
a.Size(), res,
|
||||
[=] MFEM_HOST_DEVICE(int i, unsigned &r) { r |= dptr[i]; },
|
||||
BOrReducer<unsigned> {}, use_dev, workspace);
|
||||
BOrReducer<unsigned> {}, use_dev);
|
||||
CAPTURE(use_dev);
|
||||
REQUIRE(res == 0x3ffu);
|
||||
}
|
||||
@@ -165,7 +161,6 @@ TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
|
||||
|
||||
TEST_CASE("Reduce Min", "[Reduction],[GPU]")
|
||||
{
|
||||
Array<int> workspace;
|
||||
Array<int> a(1000);
|
||||
auto hptr = a.HostReadWrite();
|
||||
for (int i = 0; i < a.Size(); ++i)
|
||||
@@ -190,7 +185,7 @@ TEST_CASE("Reduce Min", "[Reduction],[GPU]")
|
||||
r = dptr[i];
|
||||
}
|
||||
},
|
||||
MinReducer<int> {}, use_dev, workspace);
|
||||
MinReducer<int> {}, use_dev);
|
||||
CAPTURE(use_dev);
|
||||
REQUIRE(res == -10);
|
||||
}
|
||||
@@ -198,7 +193,6 @@ TEST_CASE("Reduce Min", "[Reduction],[GPU]")
|
||||
|
||||
TEST_CASE("Reduce Max", "[Reduction],[GPU]")
|
||||
{
|
||||
Array<int> workspace;
|
||||
Array<int> a(1000);
|
||||
auto hptr = a.HostReadWrite();
|
||||
for (int i = 0; i < a.Size(); ++i)
|
||||
@@ -223,7 +217,7 @@ TEST_CASE("Reduce Max", "[Reduction],[GPU]")
|
||||
r = dptr[i];
|
||||
}
|
||||
},
|
||||
MaxReducer<int> {}, use_dev, workspace);
|
||||
MaxReducer<int> {}, use_dev);
|
||||
CAPTURE(use_dev);
|
||||
REQUIRE(res == 999 - 10);
|
||||
}
|
||||
@@ -231,7 +225,6 @@ TEST_CASE("Reduce Max", "[Reduction],[GPU]")
|
||||
|
||||
TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
|
||||
{
|
||||
Array<DevicePair<int, int>> workspace;
|
||||
Array<int> a(1000);
|
||||
auto hptr = a.HostReadWrite();
|
||||
for (int i = 0; i < a.Size(); ++i)
|
||||
@@ -262,7 +255,7 @@ TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
|
||||
r.second = dptr[i];
|
||||
}
|
||||
},
|
||||
MinMaxReducer<int> {}, use_dev, workspace);
|
||||
MinMaxReducer<int> {}, use_dev);
|
||||
CAPTURE(use_dev);
|
||||
REQUIRE(res.first == -10);
|
||||
REQUIRE(res.second == a.Size() - 11);
|
||||
@@ -271,7 +264,6 @@ TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
|
||||
|
||||
TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
|
||||
{
|
||||
Array<DevicePair<double, int>> workspace;
|
||||
Array<double> a(1000);
|
||||
auto hptr = a.HostReadWrite();
|
||||
for (int i = 0; i < a.Size(); ++i)
|
||||
@@ -297,7 +289,7 @@ TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
|
||||
r.second = i;
|
||||
}
|
||||
},
|
||||
ArgMinReducer<double, int> {}, use_dev, workspace);
|
||||
ArgMinReducer<double, int> {}, use_dev);
|
||||
CAPTURE(use_dev);
|
||||
REQUIRE(res.first == -10);
|
||||
REQUIRE(res.second >= 0);
|
||||
@@ -308,7 +300,6 @@ TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
|
||||
|
||||
TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
|
||||
{
|
||||
Array<DevicePair<double, int>> workspace;
|
||||
Array<double> a(1000);
|
||||
|
||||
auto hptr = a.HostReadWrite();
|
||||
@@ -337,7 +328,7 @@ TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
|
||||
r.second = i;
|
||||
}
|
||||
},
|
||||
ArgMaxReducer<double, int> {}, use_dev, workspace);
|
||||
ArgMaxReducer<double, int> {}, use_dev);
|
||||
CAPTURE(use_dev);
|
||||
REQUIRE(res.first == a.Size() - 11);
|
||||
REQUIRE(res.second >= 0);
|
||||
@@ -348,7 +339,6 @@ TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
|
||||
|
||||
TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
|
||||
{
|
||||
Array<MinMaxLocScalar<double, int>> workspace;
|
||||
Array<double> a(1000);
|
||||
auto hptr = a.HostReadWrite();
|
||||
for (int i = 0; i < a.Size(); ++i)
|
||||
@@ -383,7 +373,7 @@ TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
|
||||
r.max_loc = i;
|
||||
}
|
||||
},
|
||||
ArgMinMaxReducer<double, int> {}, use_dev, workspace);
|
||||
ArgMinMaxReducer<double, int> {}, use_dev);
|
||||
CAPTURE(use_dev);
|
||||
REQUIRE(res.min_val == -10);
|
||||
REQUIRE(res.min_loc >= 0);
|
||||
|
||||
Reference in New Issue
Block a user