65 using namespace navier;
71 double t_final = 250*dt;
74 double fluid_kin_vis = 0.001;
77 int ode_solver_type = 3;
78 double alpha = 1.0e-2;
120 virtual void ImplicitSolve(
const double dt,
const Vector &
u,
Vector &k);
125 virtual ~ConductionOperator();
130 int x = 0,
int y = 0,
int w = 400,
int h = 400,
133 int main(
int argc,
char *argv[])
136 Mpi::Init(argc, argv);
137 int myid = Mpi::WorldRank();
144 rs_levels(lim_meshes);
147 bool visualization =
true;
150 args.
AddOption(&np_list[0],
"-np1",
"--np1",
151 "number of MPI ranks for mesh 1");
152 args.
AddOption(&np_list[1],
"-np2",
"--np2",
153 "number of MPI ranks for mesh 1");
154 args.
AddOption(&rs_levels[0],
"-r1",
"--refine-serial 1",
155 "Number of times to refine the mesh 1 uniformly in serial.");
156 args.
AddOption(&rs_levels[1],
"-r2",
"--refine-serial 2",
157 "Number of times to refine the mesh 2 uniformly in serial.");
158 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
159 "--no-visualization",
160 "Enable or disable GLVis visualization.");
173 const int nmeshes = 2;
174 mesh_file_list[0] =
"fluid-cht.mesh";
175 mesh_file_list[1] =
"solid-cht.mesh";
178 MPI_Comm *comml =
new MPI_Comm;
181 for (
int i = 0; i < nmeshes; i++)
184 if (myid < npsum) { color = i;
break; }
186 MPI_Comm_split(MPI_COMM_WORLD, color, myid, comml);
187 int myidlocal, numproclocal;
188 MPI_Comm_rank(*comml, &myidlocal);
189 MPI_Comm_size(*comml, &numproclocal);
191 Mesh *mesh =
new Mesh(mesh_file_list[color], 1, 1);
195 for (
int lev = 0; lev < rs_levels[color]; lev++)
201 if (color == 0 && myidlocal == 0)
203 std::cout <<
"Number of elements: " << mesh->
GetNE() << std::endl;
208 pmesh =
new ParMesh(*comml, *mesh);
218 ConductionOperator *coper = NULL;
224 bool last_step =
false;
246 flowsolver->
Setup(dt);
255 switch (
schwarz.ode_solver_type)
263 case 12: ode_solver =
new RK2Solver(0.5);
break;
265 case 14: ode_solver =
new RK4Solver;
break;
272 std::cout <<
"Unknown ODE solver type: " <<
schwarz.ode_solver_type <<
'\n';
295 finder.
Setup(*pmesh, color);
300 for (
int i = 0; i < color_array.
Size(); i++)
302 color_array[i] = (
unsigned int)color;
309 finder.
Interpolate(vxyz, color_array, *u_gf, interp_vals);
317 coper =
new ConductionOperator(*fes_s,
schwarz.alpha,
schwarz.kappa,
319 coper->SetParameters(adv_gf_c);
326 int Ww = 350,
Wh = 350;
327 int Wx = color*
Ww+10,
Wy = 0;
342 if (ode_solver) { ode_solver->
Init(*coper); }
344 for (
int step = 0; !last_step; ++step)
346 if (
t + dt >= t_final - dt / 2)
354 flowsolver->
Step(
t, dt, step);
359 ode_solver->
Step(t_tdof,
t, dt);
362 finder.
Interpolate(vxyz, color_array, *u_gf, interp_vals);
367 coper->SetParameters(adv_gf_c);
384 if (color == 0 && myidlocal == 0)
386 printf(
"%11s %11s %11s\n",
"Time",
"dt",
"CFL");
387 printf(
"%.5E %.5E %.5E\n",
t, dt,cfl);
397 if (color == 1) {
delete u_gf; }
413 T(NULL), current_dt(0.0),
414 M_solver(
f.GetComm()), T_solver(
f.GetComm()), udir(10), z(height)
416 const double rel_tol = 1e-8;
418 Array<int> ess_bdr(
f.GetParMesh()->bdr_attributes.Max());
425 f.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
430 M->FormSystemMatrix(ess_tdof_list, Mmat);
432 M_solver.iterative_mode =
false;
433 M_solver.SetRelTol(rel_tol);
434 M_solver.SetAbsTol(0.0);
435 M_solver.SetMaxIter(100);
436 M_solver.SetPrintLevel(0);
437 M_prec.SetType(HypreSmoother::Jacobi);
438 M_solver.SetPreconditioner(M_prec);
439 M_solver.SetOperator(Mmat);
444 T_solver.iterative_mode =
false;
445 T_solver.SetRelTol(rel_tol);
446 T_solver.SetAbsTol(0.0);
447 T_solver.SetMaxIter(100);
448 T_solver.SetPrintLevel(0);
449 T_solver.SetPreconditioner(T_prec);
451 SetParameters(adv_gf_c);
462 K->EliminateVDofsInRHS(ess_tdof_list,
u, z);
464 M_solver.Mult(z, du_dt);
469 void ConductionOperator::ImplicitSolve(
const double dt,
477 T =
Add(1.0, Mmat, dt, Kmat);
479 T_solver.SetOperator(*T);
481 MFEM_VERIFY(dt == current_dt,
"");
484 K->EliminateVDofsInRHS(ess_tdof_list,
u, z);
486 T_solver.Mult(z, du_dt);
494 u_alpha_gf.ProjectCoefficient(kapfuncoef);
504 K->FormSystemMatrix(ess_tdof_list, Kmat);
509 ConductionOperator::~ConductionOperator()
518 int x,
int y,
int w,
int h,
bool vec)
522 MPI_Comm comm = pmesh.
GetComm();
525 MPI_Comm_size(comm, &num_procs);
526 MPI_Comm_rank(comm, &myid);
528 bool newly_opened =
false;
529 int connection_failed;
541 sock <<
"solution\n";
547 if (myid == 0 && newly_opened)
550 ?
"mAcRjlmm" :
"mmaaAcl";
552 sock <<
"window_title '" << title <<
"'\n" 553 <<
"window_geometry " 554 << x <<
" " << y <<
" " << w <<
" " << h <<
"\n" 556 if ( vec ) { sock <<
"vvv"; }
562 connection_failed = !sock && !newly_opened;
564 MPI_Bcast(&connection_failed, 1, MPI_INT, 0, comm);
566 while (connection_failed);
578 if (std::fabs(xi+2.5)<1.e-5) {
u(0) = 0.25*yi*(3-yi)/(1.5*1.5); }
585 return x(1) <= 1.0 && std::fabs(x(0)) < 0.5 ? 5.: 1.0;
594 t_init = 10*(std::exp(-x(1)*x(1)));
596 if (std::fabs(x(0)) >= 0.5)
598 double dx = std::fabs(x(0))-0.5;
599 t_init *= std::exp(-10*dx*dx);
void PrintTimingData()
Print timing summary of the solving routine.
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Conjugate gradient method.
double kappa_fun(const Vector &x)
Solid data.
void PrintOptions(std::ostream &out) const
Print the options.
int Dimension() const
Dimension of the reference space used within the elements.
void SaveAsOne(const char *fname, int precision=16) const
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
void PrintUsage(std::ostream &out) const
Print the usage message.
Base abstract class for first order time dependent operators.
Vector coefficient that is constant in space and time.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
virtual void Step(Vector &x, double &t, double &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
int Size() const
Returns the size of the vector.
void vel_dbc(const Vector &x, double t, Vector &u)
Fluid data.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
void Interpolate(const Vector &point_pos, Array< unsigned int > &point_id, const GridFunction &field_in, Vector &field_out, int point_pos_ordering=Ordering::byNODES)
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
bool Good() const
Return true if the command line options were parsed successfully.
Abstract parallel finite element space.
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Backward Euler ODE solver. L-stable.
struct schwarz_common schwarz
void Setup(double dt)
Initialize forms, solvers and preconditioners.
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
double ComputeCFL(ParGridFunction &u, double dt)
Compute CFL.
OversetFindPointsGSLIB enables use of findpts for arbitrary number of overlapping grids...
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
ParMesh * GetParMesh() const
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
void EnablePA(bool pa)
Enable partial assembly for every operator.
void AddVelDirichletBC(VectorCoefficient *coeff, Array< int > &attr)
Add a Dirichlet boundary condition to the velocity field.
Parallel smoothers in hypre.
void VisualizeField(socketstream &sock, const char *vishost, int visport, ParGridFunction &gf, const char *title, int x=0, int y=0, int w=400, int h=400, bool vec=false)
FiniteElementSpace * FESpace()
int main(int argc, char *argv[])
Mesh * GetMesh() const
Returns the mesh.
The classical explicit forth-order Runge-Kutta method, RK4.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
bool is_open()
True if the socketstream is open, false otherwise.
ParFiniteElementSpace * ParFESpace() const
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set 'var' to receive the value. Enable/disable tags are used to set the bool...
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Third-order, strong stability preserving (SSP) Runge-Kutta method.
void PrintAsOne(std::ostream &out=mfem::out) const
ParGridFunction * GetCurrentVelocity()
Return a pointer to the current velocity ParGridFunction.
int GetNE() const
Returns number of elements.
Implicit midpoint method. A-stable, not L-stable.
void SetGridFunction(const GridFunction *gf)
Set the grid function for this coefficient. Also sets the Vector dimension to match that of the gf...
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
int Size() const
Return the logical size of the array.
double temp_init(const Vector &x)
A general function coefficient.
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Vector coefficient defined by a vector GridFunction.
Arbitrary order H1-conforming (continuous) finite elements.
void GetNodes(Vector &node_coord) const
void Step(double &time, double dt, int cur_step, bool provisional=false)
Compute solution at the next time step t+dt.
double u(const Vector &xvec)
Class for parallel grid function.
The classical forward Euler method.
Wrapper for hypre's ParCSR matrix class.
Class for parallel meshes.
Transient incompressible Navier Stokes solver in a split scheme formulation.
void Setup(Mesh &m, const int meshid, GridFunction *gfmax=NULL, const double bb_t=0.1, const double newt_tol=1.0e-12, const int npt_max=256)
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
double f(const Vector &p)