42 int main(
int argc,
char *argv[])
48 MPI_Init(&argc, &argv);
49 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
50 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
51 bool verbose = (myid == 0);
54 const char *mesh_file =
"../data/star.mesh";
56 bool visualization = 1;
59 args.
AddOption(&mesh_file,
"-m",
"--mesh",
62 "Finite element order (polynomial degree).");
63 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
65 "Enable or disable GLVis visualization.");
85 ifstream imesh(mesh_file);
90 cerr <<
"\nCan not open mesh file: " << mesh_file <<
'\n' << endl;
95 mesh =
new Mesh(imesh, 1, 1);
105 (int)floor(log(10000./mesh->
GetNE())/log(2.)/
dim);
106 for (
int l = 0; l < ref_levels; l++)
118 int par_ref_levels = 2;
119 for (
int l = 0; l < par_ref_levels; l++)
138 std::cout <<
"***********************************************************\n";
139 std::cout <<
"dim(R) = " << dimR <<
"\n";
140 std::cout <<
"dim(W) = " << dimW <<
"\n";
141 std::cout <<
"dim(R+W) = " << dimR + dimW <<
"\n";
142 std::cout <<
"***********************************************************\n";
151 block_offsets[0] = 0;
152 block_offsets[1] = R_space->
GetVSize();
153 block_offsets[2] = W_space->
GetVSize();
157 block_trueOffsets[0] = 0;
158 block_trueOffsets[1] = R_space->
TrueVSize();
159 block_trueOffsets[2] = W_space->
TrueVSize();
175 BlockVector trueX(block_trueOffsets), trueRhs(block_trueOffsets);
265 solver.
Mult(trueRhs, trueX);
272 <<
" iterations with a residual norm of " << solver.
GetFinalNorm() <<
".\n";
275 <<
" iterations. Residual norm is " << solver.
GetFinalNorm() <<
".\n";
276 std::cout <<
"MINRES solver took " << chrono.
RealTime() <<
"s. \n";
284 u->
Update(R_space, x.GetBlock(0), 0);
285 p->
Update(W_space, x.GetBlock(1), 0);
289 int order_quad = max(2, 2*order+1);
291 for (
int i=0; i < Geometry::NumGeom; ++i)
303 std::cout <<
"|| u_h - u_ex || / || u_ex || = " << err_u / norm_u <<
"\n";
304 std::cout <<
"|| p_h - p_ex || / || p_ex || = " << err_p / norm_p <<
"\n";
310 ostringstream mesh_name, u_name, p_name;
311 mesh_name <<
"mesh." << setfill(
'0') << setw(6) << myid;
312 u_name <<
"sol_u." << setfill(
'0') << setw(6) << myid;
313 p_name <<
"sol_p." << setfill(
'0') << setw(6) << myid;
315 ofstream mesh_ofs(mesh_name.str().c_str());
316 mesh_ofs.precision(8);
317 pmesh->
Print(mesh_ofs);
319 ofstream u_ofs(u_name.str().c_str());
323 ofstream p_ofs(p_name.str().c_str());
337 char vishost[] =
"localhost";
340 u_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n";
342 u_sock <<
"solution\n" << *pmesh << *u <<
"window_title 'Velocity'"
348 p_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n";
350 p_sock <<
"solution\n" << *pmesh << *p <<
"window_title 'Pressure'"
393 u(0) = - exp(xi)*sin(yi)*cos(zi);
394 u(1) = - exp(xi)*cos(yi)*cos(zi);
398 u(2) = exp(xi)*sin(yi)*sin(zi);
414 return exp(xi)*sin(yi)*cos(zi);
Class for domain integration L(v) := (f, v)
Class for integration rule.
int GetNumIterations() const
double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
void fFun(const Vector &x, Vector &f)
virtual void RegisterField(const char *field_name, GridFunction *gf)
Add a grid function to the collection and update the root file.
Subclass constant coefficient.
HYPRE_Int GlobalTrueVSize()
int Size() const
Returns the size of the vector.
int GetNE() const
Returns number of elements.
virtual void Save()
Save the collection and a VisIt root file.
virtual void Save(std::ostream &out) const
Abstract parallel finite element space.
void SetBlock(int iRow, int iCol, Operator *op)
Add a block op in the block-entry (iblock, jblock).
bool iterative_mode
If true, use the second argument of Mult as an initial guess.
double GetFinalNorm() const
The BoomerAMG solver in hypre.
virtual void Mult(const Vector &b, Vector &x) const
Operator application.
HYPRE_Int GetGlobalNumRows() const
IntegrationRules IntRules(0)
A global object with all integration rules (defined in intrules.cpp)
void SetPrintLevel(int print_lvl)
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
Jacobi preconditioner in hypre.
double f_natural(const Vector &x)
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Data collection with VisIt I/O routines.
void SetMaxIter(int max_it)
double pFun_ex(const Vector &x)
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
void PrintUsage(std::ostream &out) const
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Wrapper for hypre's parallel vector class.
void SetAbsTol(double atol)
int main(int argc, char *argv[])
void SetRelTol(double rtol)
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
void PartialSum()
Partial Sum.
double ComputeGlobalLpNorm(double p, Coefficient &coeff, ParMesh &pmesh, const IntegrationRule *irs[])
void Distribute(const Vector *tv)
HypreParMatrix * Transpose()
Returns the transpose of *this.
void uFun_ex(const Vector &x, Vector &u)
double gFun(const Vector &x)
void InvScaleRows(const Vector &s)
Scale the local row i by 1./s(i)
void PrintOptions(std::ostream &out) const
Abstract class for hypre's solvers and preconditioners.
for VectorFiniteElements (Nedelec, Raviart-Thomas)
class for C-function coefficient
HYPRE_Int * GetRowStarts() const
Class for parallel grid function.
Wrapper for hypre's ParCSR matrix class.
A class to handle Block systems in a matrix-free implementation.
Vector & GetBlock(int i)
Get the i-th vector in the block.
Class for parallel meshes.
Integrator for (Q u, v) for VectorFiniteElements.
HypreParMatrix * ParMult(HypreParMatrix *A, HypreParMatrix *B)
Returns the matrix A * B.
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
Arbitrary order "L2-conforming" discontinuous finite elements.
virtual void Print(std::ostream &out=std::cout) const