Commit a40c59bb authored by grandper's avatar grandper
Browse files

Update of the init policies

parent 3ec248e7
......@@ -88,7 +88,8 @@ using namespace LifeV;
namespace
{
typedef RegionMesh<LinearTetra> mesh_Type;
//typedef RegionMesh<LinearTetra> mesh_Type;
typedef RegionMesh<QuadraticTetra> mesh_Type;
typedef MatrixEpetra<Real> matrix_Type;
typedef VectorEpetra vector_Type;
typedef boost::shared_ptr<VectorEpetra> vectorPtr_Type;
......@@ -104,9 +105,9 @@ typedef TimeIterationPolicyLinear< AssemblyPolicyNavierStokesSemiImplicit > Semi
typedef TimeIterationPolicyNonlinearIncremental< AssemblyPolicyNavierStokesNewton > Newton;
typedef TimeIterationPolicyNonlinearIncremental< AssemblyPolicyNavierStokesPicard > Picard;
typedef TimeIterationPolicyNonlinear< AssemblyPolicyNavierStokesPicard > PicardOseen;
typedef InitPolicySolver<Stokes> InitStokes;
typedef InitPolicySolver<GStokes> InitGStokes;
typedef InitPolicyInterpolation InitInter;
typedef InitPolicySolver< mesh_Type, Stokes > InitStokes;
typedef InitPolicySolver< mesh_Type, GStokes > InitGStokes;
typedef InitPolicyInterpolation< mesh_Type > InitInter;
typedef InitPolicyProjection<SolverPolicyLinearSolver> InitProj;
typedef ExporterPolicyNoExporter NoExporter;
typedef ExporterPolicyHDF5 HDF5Exporter;
......@@ -200,7 +201,7 @@ main ( int argc, char** argv )
std::string meshPath = problemList.get ( "Resources path", "./Resources" );
meshPath.append ("/");
boost::shared_ptr<NavierStokesProblem<mesh_Type> > nsProblem;
boost::shared_ptr< NavierStokesProblem< mesh_Type > > nsProblem;
if ( benchmark == "Cavity" )
{
nsProblem.reset ( new NavierStokesCavity );
......
......@@ -24,8 +24,6 @@ SET(NavierStokesSolver_SOURCES
solver/NavierStokesSolver/AssemblyPolicyNavierStokesNewton.cpp
solver/NavierStokesSolver/AssemblyPolicyNavierStokesPicard.cpp
solver/NavierStokesSolver/ExporterPolicyHDF5.cpp
solver/NavierStokesSolver/InitPolicyInterpolation.cpp
solver/NavierStokesSolver/NavierStokesProblem.cpp
solver/NavierStokesSolver/SolverPolicyLinearSolver.cpp
CACHE INTERNAL "")
......
//@HEADER
/*
*******************************************************************************
Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
This file is part of LifeV.
LifeV is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
LifeV is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with LifeV. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************
*/
//@HEADER
/*!
@file InitPolicyInterpolation class
@brief This class is a strategy to initialize a Navier-Stokes problem
@author Gwenol Grandperrin <gwenol.grandperrin@epfl.ch>
@date 11-12-2012
*/
#include <lifev/navier_stokes/solver/NavierStokesSolver/InitPolicyInterpolation.hpp>
#include <string>
#include <lifev/core/util/LifeChrono.hpp>
namespace LifeV
{
void
InitPolicyInterpolation::
setupInit ( Teuchos::ParameterList& /*list*/ )
{
}
void
InitPolicyInterpolation::
initSimulation ( bcContainerPtr_Type /*bchandler*/,
vectorPtr_Type solution )
{
ASSERT ( problem()->hasExactSolution(), "Error: You cannot use the interpolation method if the problem has not an analytical solution." );
Real currentTime = initialTime() - timestep() * bdf()->order();
UInt pressureOffset = uFESpace()->fieldDim() * uFESpace()->dof().numTotalDof();
vectorPtr_Type velocity;
velocity.reset ( new vector_Type ( uFESpace()->map(), Unique ) );
vectorPtr_Type pressure;
pressure.reset ( new vector_Type ( pFESpace()->map(), Unique ) );
*solution = 0.0;
uFESpace()->interpolate ( problem()->uexact(), *velocity, currentTime );
pFESpace()->interpolate ( problem()->pexact(), *pressure, currentTime );
solution->add ( *velocity );
solution->add ( *pressure, pressureOffset );
bdf()->setInitialCondition ( *solution );
currentTime += timestep();
for ( ; currentTime <= initialTime() + timestep() / 2.0; currentTime += timestep() )
{
*solution = 0.0;
*velocity = 0.0;
*pressure = 0.0;
uFESpace()->interpolate ( problem()->uexact(), *velocity, currentTime );
pFESpace()->interpolate ( problem()->pexact(), *pressure, currentTime );
solution->add ( *velocity );
solution->add ( *pressure, pressureOffset );
// Updating bdf
bdf()->shiftRight ( *solution );
}
}
} // namespace LifeV
......@@ -36,6 +36,7 @@
#define INITPOLICYINTERPOLATION_HPP
#include <iostream>
#include <string>
#include <boost/shared_ptr.hpp>
// Tell the compiler to ignore specific kind of warnings:
......@@ -57,6 +58,7 @@
#pragma GCC diagnostic warning "-Wunused-parameter"
#include <lifev/core/LifeV.hpp>
#include <lifev/core/util/LifeChrono.hpp>
#include <lifev/core/array/VectorEpetra.hpp>
#include <lifev/core/mesh/RegionMesh.hpp>
#include <lifev/core/fem/FESpace.hpp>
......@@ -67,13 +69,13 @@
namespace LifeV
{
template< class mesh_Type >
struct InitPolicyInterpolation
{
typedef VectorEpetra vector_Type;
typedef boost::shared_ptr<VectorEpetra> vectorPtr_Type;
typedef MapEpetra map_Type;
typedef boost::shared_ptr<map_Type> mapPtr_Type;
typedef RegionMesh<LinearTetra> mesh_Type;
typedef FESpace< mesh_Type, map_Type > fespace_Type;
typedef boost::shared_ptr< fespace_Type > fespacePtr_Type;
typedef TimeAdvanceBDF<vector_Type> bdf_Type;
......@@ -95,6 +97,55 @@ struct InitPolicyInterpolation
virtual NSProblemPtr_Type problem() const = 0;
};
template< class mesh_Type >
void
InitPolicyInterpolation< mesh_Type >::
setupInit ( Teuchos::ParameterList& /*list*/ )
{
}
template< class mesh_Type >
void
InitPolicyInterpolation< mesh_Type >::
initSimulation ( bcContainerPtr_Type /*bchandler*/,
vectorPtr_Type solution )
{
ASSERT ( problem()->hasExactSolution(), "Error: You cannot use the interpolation method if the problem has not an analytical solution." );
Real currentTime = initialTime() - timestep() * bdf()->order();
UInt pressureOffset = uFESpace()->fieldDim() * uFESpace()->dof().numTotalDof();
vectorPtr_Type velocity;
velocity.reset ( new vector_Type ( uFESpace()->map(), Unique ) );
vectorPtr_Type pressure;
pressure.reset ( new vector_Type ( pFESpace()->map(), Unique ) );
*solution = 0.0;
uFESpace()->interpolate ( problem()->uexact(), *velocity, currentTime );
pFESpace()->interpolate ( problem()->pexact(), *pressure, currentTime );
solution->add ( *velocity );
solution->add ( *pressure, pressureOffset );
bdf()->setInitialCondition ( *solution );
currentTime += timestep();
for ( ; currentTime <= initialTime() + timestep() / 2.0; currentTime += timestep() )
{
*solution = 0.0;
*velocity = 0.0;
*pressure = 0.0;
uFESpace()->interpolate ( problem()->uexact(), *velocity, currentTime );
pFESpace()->interpolate ( problem()->pexact(), *pressure, currentTime );
solution->add ( *velocity );
solution->add ( *pressure, pressureOffset );
// Updating bdf
bdf()->shiftRight ( *solution );
}
}
} // namespace LifeV
#endif /* INITPOLICYINTERPOLATION_HPP */
......@@ -73,14 +73,13 @@
namespace LifeV
{
template< class SolverPolicy = SolverPolicyLinearSolver >
template< class mesh_Type, class SolverPolicy = SolverPolicyLinearSolver >
struct InitPolicyProjection : public virtual SolverPolicy, public AssemblyPolicyStokes
{
typedef VectorEpetra vector_Type;
typedef boost::shared_ptr<vector_Type> vectorPtr_Type;
typedef MatrixEpetra<Real> matrix_Type;
typedef boost::shared_ptr<matrix_Type> matrixPtr_Type;
typedef RegionMesh<LinearTetra> mesh_Type;
typedef MeshPartitioner< mesh_Type > meshPartitioner_Type;
typedef MapEpetra map_Type;
typedef boost::shared_ptr<map_Type> mapPtr_Type;
......@@ -113,18 +112,18 @@ private:
};
template< class SolverPolicy >
template< class mesh_Type, class SolverPolicy >
void
InitPolicyProjection<SolverPolicy>::
InitPolicyProjection< mesh_Type, SolverPolicy >::
setupInit ( Teuchos::ParameterList& list )
{
Teuchos::ParameterList assemblyList = list.sublist ( "Assembly: Parameter list" );
M_list = list;
}
template< class SolverPolicy >
template< class mesh_Type, class SolverPolicy >
void
InitPolicyProjection<SolverPolicy>::
InitPolicyProjection< mesh_Type, SolverPolicy >::
initSimulation ( bcContainerPtr_Type bchandler,
vectorPtr_Type solution )
{
......
......@@ -68,12 +68,11 @@
namespace LifeV
{
template< class TimeIterationPolicy >
template< class mesh_Type, class TimeIterationPolicy >
struct InitPolicySolver : public virtual TimeIterationPolicy
{
typedef VectorEpetra vector_Type;
typedef boost::shared_ptr<VectorEpetra> vectorPtr_Type;
typedef RegionMesh<LinearTetra> mesh_Type;
typedef MeshPartitioner< mesh_Type > meshPartitioner_Type;
typedef MapEpetra map_Type;
typedef boost::shared_ptr<map_Type> mapPtr_Type;
......@@ -97,17 +96,17 @@ struct InitPolicySolver : public virtual TimeIterationPolicy
virtual bdfPtr_Type bdf() const = 0;
};
template< class TimeIterationPolicy >
template< class mesh_Type, class TimeIterationPolicy >
void
InitPolicySolver< TimeIterationPolicy >::
InitPolicySolver< mesh_Type, TimeIterationPolicy >::
setupInit ( Teuchos::ParameterList& list )
{
TimeIterationPolicy::initTimeIteration ( list );
}
template< class TimeIterationPolicy >
template< class mesh_Type, class TimeIterationPolicy >
void
InitPolicySolver< TimeIterationPolicy >::
InitPolicySolver< mesh_Type, TimeIterationPolicy >::
initSimulation ( bcContainerPtr_Type bchandler,
vectorPtr_Type solution )
{
......
......@@ -116,7 +116,7 @@ public:
//@{
//! Getter for the problem mesh
virtual void mesh ( boost::shared_ptr<mesh_Type>& mesh ) const = 0;
virtual void mesh ( boost::shared_ptr< mesh_Type >& mesh ) const = 0;
//! Getter for the boundary conditions in the provided BCHandler
/*!
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment