40 #include <Epetra_ConfigDefs.h> 43 #include <Epetra_MpiComm.h> 45 #include <Epetra_SerialComm.h> 49 #include <lifev/core/array/MapEpetra.hpp> 50 #include <lifev/core/mesh/MeshData.hpp> 51 #include <lifev/core/mesh/MeshPartitioner.hpp> 52 #include <lifev/core/algorithm/SolverAztecOO.hpp> 53 #include <lifev/core/algorithm/LinearSolver.hpp> 54 #include <lifev/core/fem/Assembly.hpp> 55 #include <lifev/core/fem/AssemblyElemental.hpp> 56 #include <lifev/core/fem/BCManage.hpp> 57 #include <lifev/core/fem/FESpace.hpp> 58 #include <lifev/core/fem/TimeAdvanceBDFVariableStep.hpp> 60 #include <lifev/core/algorithm/PreconditionerIfpack.hpp> 61 #include <lifev/core/algorithm/PreconditionerML.hpp> 63 #include <lifev/core/filter/Exporter.hpp> 64 #include <lifev/core/filter/ExporterEmpty.hpp> 65 #include <lifev/core/filter/ExporterEnsight.hpp> 66 #include <lifev/core/filter/ExporterHDF5.hpp> 67 #include <lifev/core/mesh/RegionMesh.hpp> 69 #include <Teuchos_ParameterList.hpp> 70 #include <Teuchos_XMLParameterListHelpers.hpp> 71 #include <Teuchos_RCP.hpp> 79 using namespace LifeV;
111 test_bdf::test_bdf (
int argc,
char** argv) :
112 Members (
new Private)
115 const std::string data_file_name =
116 command_line.follow (
"data_2d", 2,
"-f",
"--file");
117 GetPot dataFile (data_file_name);
119 Members->data_file_name = data_file_name;
122 std::cout <<
"Epetra Initialization" << std::endl;
123 Members->comm.reset (
new Epetra_MpiComm (MPI_COMM_WORLD) );
125 Members->comm.reset (
new Epetra_SerialComm() );
135 typedef VectorEpetra vector_Type;
136 typedef std::shared_ptr<vector_Type> vectorPtr_Type;
138 typedef LifeV::Preconditioner basePrec_Type;
139 typedef std::shared_ptr<basePrec_Type> basePrecPtr_Type;
140 typedef LifeV::PreconditionerIfpack prec_Type;
141 typedef std::shared_ptr<prec_Type> precPtr_Type;
144 GetPot dataFile (Members->data_file_name.c_str() );
146 bool verbose = (Members->comm->MyPID() == 0);
150 std::cout <<
"The BDF Solver" << std::flush;
156 BCFunctionBase g_Ess (AnalyticalSol_2d::u);
160 bc.addBC (
"Top", TOP, Essential, Full, g_Ess, 1);
161 bc.addBC (
"Bottom", BOTTOM, Essential, Full, g_Ess, 1);
162 bc.addBC (
"Left", LEFT, Essential, Full, g_Ess, 1);
163 bc.addBC (
"Right", RIGHT, Essential, Full, g_Ess, 1);
166 Members->comm->Barrier();
167 MeshData meshData (dataFile, (
"bdf/" + discretization_section).c_str() );
168 std::shared_ptr<regionMesh> fullMeshPtr (
new regionMesh ( Members->comm ) );
169 readMesh (*fullMeshPtr, meshData);
170 std::shared_ptr<regionMesh> meshPtr;
173 meshPtr = meshPart.meshPartition();
178 std::shared_ptr<FESpace<regionMesh, MapEpetra> > feSpacePtr (
179 new FESpace<regionMesh, MapEpetra> (
180 meshPtr, dataFile ( (
"bdf/" + discretization_section +
"/order").c_str(),
"P2"), 1, Members->comm) );
183 std::cout <<
" Number of unknowns : " 184 << feSpacePtr->map().map (Unique)->NumGlobalElements()
187 bc.bcUpdate (* (feSpacePtr->mesh() ), feSpacePtr->feBd(), feSpacePtr->dof() );
191 MatrixElemental elmat (feSpacePtr->fe().nbFEDof(), 1, 1);
192 MatrixEpetra<
double> matM (feSpacePtr->map() );
193 std::shared_ptr<MatrixEpetra<
double> > matA_ptr (
194 new MatrixEpetra<
double> (feSpacePtr->map() ) );
195 std::shared_ptr<VectorEpetra> u (
new VectorEpetra (feSpacePtr->map(), Unique) );
196 std::shared_ptr<VectorEpetra> f (
new VectorEpetra (feSpacePtr->map(), Unique) );
200 Members->comm->Barrier();
202 for (UInt iVol = 0; iVol < feSpacePtr->mesh()->numElements(); iVol++)
204 feSpacePtr->fe().updateJac (feSpacePtr->mesh()->element (iVol) );
206 AssemblyElemental::mass (1., elmat, feSpacePtr->fe(), 0, 0);
207 assembleMatrix (matM, elmat, feSpacePtr->fe(), feSpacePtr->fe(), feSpacePtr->dof(),
208 feSpacePtr->dof(), 0, 0, 0, 0);
210 matM.globalAssemble();
211 Members->comm->Barrier();
214 std::cout <<
"\n \n -- Mass matrix assembling time = " << chrono.diff() << std::endl
220 Real Tfin = dataFile
("bdf/endtime", 10.0
);
221 Real delta_t = dataFile
("bdf/timestep", 0.5
);
223 UInt ord_bdf = dataFile
("bdf/order", 3
);
228 bdf.setInitialCondition<Real (*) (Real, Real, Real, Real, UInt), FESpace<regionMesh, MapEpetra> > (AnalyticalSol_2d::u, *u, *feSpacePtr, t0, delta_t);
234 Members->comm->Barrier();
238 std::shared_ptr<Exporter<regionMesh> > exporter;
239 std::string
const exporterType = dataFile (
"exporter/type",
"hdf5");
242 if (exporterType.compare (
"hdf5") == 0)
244 exporter.reset (
new ExporterHDF5<regionMesh > ( dataFile,
"bdf_test" ) );
249 if (exporterType.compare (
"none") == 0)
251 exporter.reset (
new ExporterEmpty<regionMesh > ( dataFile, meshPtr,
"bdf_test", Members->comm->MyPID() ) );
255 exporter.reset (
new ExporterEnsight<regionMesh > ( dataFile, meshPtr,
"bdf_test", Members->comm->MyPID() ) );
259 exporter->setPostDir (
"./" );
260 exporter->setMeshProcId ( meshPtr, Members->comm->MyPID() );
262 std::shared_ptr<VectorEpetra> u_display_ptr (
new VectorEpetra (
263 feSpacePtr->map(), exporter->mapType() ) );
264 exporter->addVariable (ExporterData<regionMesh >::ScalarField,
"u", feSpacePtr,
265 u_display_ptr, UInt (0) );
267 exporter->postProcess (0);
272 LinearSolver linearSolver;
274 Teuchos::RCP< Teuchos::ParameterList > belosList = Teuchos::rcp (
new Teuchos::ParameterList );
275 belosList = Teuchos::getParametersFromXmlFile (
"SolverParamList.xml" );
277 linearSolver.setCommunicator ( Members->comm );
278 linearSolver.setParameters ( *belosList );
280 prec_Type* precRawPtr;
281 basePrecPtr_Type precPtr;
282 precRawPtr =
new prec_Type;
283 precRawPtr->setDataFromGetPot ( dataFile,
"bdf/prec" );
284 precPtr.reset ( precRawPtr );
285 linearSolver.setPreconditioner ( precPtr );
291 matA_ptr.reset (
new MatrixEpetra<
double> (feSpacePtr->map() ) );
293 for (Real t = t0 + delta_t; t <= Tfin; t += delta_t)
295 Members->comm->Barrier();
298 std::cout <<
"Now we are at time " << t << std::endl;
303 Real coeff = bdf.coefficientDerivative (0) / delta_t;
306 for (UInt i = 0; i < feSpacePtr->mesh()->numElements(); i++)
308 feSpacePtr->fe().updateFirstDerivQuadPt (feSpacePtr->mesh()->element (i) );
310 AssemblyElemental::mass (coeff + s, elmat, feSpacePtr->fe() );
311 AssemblyElemental::stiff (visc, elmat, feSpacePtr->fe() );
312 assembleMatrix (*matA_ptr, elmat, feSpacePtr->fe(), feSpacePtr->fe(),
313 feSpacePtr->dof(), feSpacePtr->dof(), 0, 0, 0, 0);
319 std::cout <<
"A has been constructed in " << chrono.diff() <<
"s." << std::endl;
323 *f = (matM * bdf.rhsContribution() );
324 feSpacePtr->l2ScalarProduct (sf, *f, t);
325 Members->comm->Barrier();
330 std::cout <<
"*** BC Management: " << std::endl;
334 bcManage (*matA_ptr, *f, *feSpacePtr->mesh(), feSpacePtr->dof(), bc, feSpacePtr->feBd(), tgv, t);
335 matA_ptr->globalAssemble();
339 std::cout << chrono.diff() <<
"s." << std::endl;
343 Members->comm->Barrier();
346 linearSolver.setOperator ( matA_ptr );
347 linearSolver.setReusePreconditioner (
false);
348 linearSolver.setRightHandSide ( f );
349 linearSolver.solve ( u );
354 std::cout <<
"*** Solution computed in " << chrono.diff() <<
"s." 357 Members->comm->Barrier();
360 vector_Type uComputed (*u, Repeated);
362 Real L2_Error, L2_RelError;
364 L2_Error = feSpacePtr->l2Error (AnalyticalSol_2d::u, uComputed, t, &L2_RelError);
367 std::cout <<
"Error Norm L2: " << L2_Error
368 <<
"\nRelative Error Norm L2: " << L2_RelError << std::endl;
370 Members->errorNorm = L2_Error;
377 exporter->postProcess (t);
384 bool test_bdf::check()
387 GetPot dataFile (Members->data_file_name.c_str() );
388 return Members->errorNorm < dataFile (
"errorNorms/l2Error", -10e10);
double operator()(const char *VarName, const double &Default) const
const std::string discretization_section
GetPot(const int argc_, char **argv_, const char *FieldSeparator=0x0)
std::shared_ptr< Epetra_Comm > comm
TimeAdvanceBDFVariableStep - Backward differencing formula time discretization with non-uniform time ...
std::shared_ptr< std::vector< Int > > M_isOnProc
MeshData - class for handling spatial discretization.
double Real
Generic real data.
std::string data_file_name
int operator()(const char *VarName, int Default) const
std::function< Real(Real const &, Real const &, Real const &, Real const &, ID const &) > fct_type
RegionMesh< LinearTriangle > regionMesh
uint32_type UInt
generic unsigned integer (used mainly for addressing)