Commit 37569a73 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#820 Poromechanics: added transient source in implicit step fluid. Currently...

#820 Poromechanics: added transient source in implicit step fluid. Currently fails at runtime because porosity parameter is not the same as in explicit: we need it defined on the new felt space as well...
parent a6c96fbd
......@@ -333,6 +333,43 @@ Domain11 = {
} -- Domain11
-- Inlet border
Domain12 = {
-- Index of the geometric mesh upon which the domain is defined (as defined in the present file). Might be
-- left empty if domain not limited to one mesh; at most one value is expected here.
-- Expected format: {VALUE1, VALUE2, ...}
mesh_index = { 1 },
-- List of dimensions encompassed by the domain. Might be left empty if no restriction at all upon
-- dimensions.
-- Expected format: {VALUE1, VALUE2, ...}
-- Constraint: ops_in(v, {0, 1, 2, 3})
dimension_list = { 1 },
-- List of mesh labels encompassed by the domain. Might be left empty if no restriction at all upon mesh
-- labels. This parameter does not make sense if no mesh is defined for the domain.
-- Expected format: {VALUE1, VALUE2, ...}
mesh_label_list = { 4 },
-- List of geometric element types considered in the domain. Might be left empty if no restriction upon
-- these. No constraint is applied at Ops level, as some geometric element types could be added after
-- generation of current input parameter file. Current list is below; if an incorrect value is put there it
-- will be detected a bit later when the domain object is built.
-- The known types when this file was generated are:
-- . Point1
-- . Segment2, Segment3
-- . Triangle3, Triangle6
-- . Quadrangle4, Quadrangle8, Quadrangle9
-- . Tetrahedron4, Tetrahedron10
-- . Hexahedron8, Hexahedron20, Hexahedron27.
-- Expected format: {"VALUE1", "VALUE2", ...}
geometric_element_type_list = {}
} -- Domain12
Domain20 = {
-- Index of the geometric mesh upon which the domain is defined (as defined in the present file). Might be
......@@ -465,6 +502,33 @@ FiniteElementSpace10 = {
-- FiniteElementSpace11 -- Inlet border
FiniteElementSpace11 = {
-- Index of the god of dof into which the finite element space is defined.
-- Expected format: VALUE
god_of_dof_index = 1,
-- Index of the domain onto which the finite element space is defined. This domain must be unidimensional.
-- Expected format: VALUE
domain_index = 12,
-- List of all unknowns defined in the finite element space. Unknowns here must be defined in this file as
-- an 'Unknown' block; expected name/identifier is the name given there.
-- Expected format: {"VALUE1", "VALUE2", ...}
unknown_list = { 'fluid_velocity'},
-- List of the shape function to use for each unknown;
-- Expected format: {"VALUE1", "VALUE2", ...}
-- shape_function_list = { 'P2', 'P1' },
shape_function_list = { 'P1' },
-- List of the numbering subset to use for each unknown;
-- Expected format: {VALUE1, VALUE2, ...}
numbering_subset_list = { 10 }
}
FiniteElementSpace20 = {
......
......@@ -446,6 +446,8 @@
BE50E9021C903BDF0038D16E /* None.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE50E9001C903BDF0038D16E /* None.cpp */; };
BE50E9031C903BDF0038D16E /* Viscoelasticity.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE50E9011C903BDF0038D16E /* Viscoelasticity.cpp */; };
BE51CDB41CA3FE9400B66D2A /* ImplicitStepFluidVariationalFormulation.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE51CDB11CA3FE9400B66D2A /* ImplicitStepFluidVariationalFormulation.cpp */; };
BE51CDE01CA4560100B66D2A /* TransientSource.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE51CDDD1CA4560100B66D2A /* TransientSource.cpp */; };
BE51CDE51CA4560F00B66D2A /* TransientSource.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE51CDE21CA4560F00B66D2A /* TransientSource.cpp */; };
BE5389FD1C897FE400D80749 /* Mpi.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BE5389FB1C897FE400D80749 /* Mpi.hpp */; };
BE5389FE1C897FE400D80749 /* Mpi.hxx in Headers */ = {isa = PBXBuildFile; fileRef = BE5389FC1C897FE400D80749 /* Mpi.hxx */; };
BE538A041C89B97800D80749 /* CoordIndexes.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE538A011C89B97800D80749 /* CoordIndexes.cpp */; };
......@@ -4415,6 +4417,12 @@
BE51CDB11CA3FE9400B66D2A /* ImplicitStepFluidVariationalFormulation.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = ImplicitStepFluidVariationalFormulation.cpp; sourceTree = "<group>"; };
BE51CDB21CA3FE9400B66D2A /* ImplicitStepFluidVariationalFormulation.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = ImplicitStepFluidVariationalFormulation.hpp; sourceTree = "<group>"; };
BE51CDB31CA3FE9400B66D2A /* ImplicitStepFluidVariationalFormulation.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = ImplicitStepFluidVariationalFormulation.hxx; sourceTree = "<group>"; };
BE51CDDD1CA4560100B66D2A /* TransientSource.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = TransientSource.cpp; sourceTree = "<group>"; };
BE51CDDE1CA4560100B66D2A /* TransientSource.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = TransientSource.hpp; sourceTree = "<group>"; };
BE51CDDF1CA4560100B66D2A /* TransientSource.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = TransientSource.hxx; sourceTree = "<group>"; };
BE51CDE21CA4560F00B66D2A /* TransientSource.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = TransientSource.cpp; sourceTree = "<group>"; };
BE51CDE31CA4560F00B66D2A /* TransientSource.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = TransientSource.hpp; sourceTree = "<group>"; };
BE51CDE41CA4560F00B66D2A /* TransientSource.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = TransientSource.hxx; sourceTree = "<group>"; };
BE52B1A71C8DDC9A009B7857 /* ModelInitialize.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = ModelInitialize.hxx; sourceTree = "<group>"; };
BE5389FB1C897FE400D80749 /* Mpi.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = Mpi.hpp; sourceTree = "<group>"; };
BE5389FC1C897FE400D80749 /* Mpi.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = Mpi.hxx; sourceTree = "<group>"; };
......@@ -7005,6 +7013,9 @@
BE60DBFB1C7F417C007B334C /* Porosity.cpp */,
BE60DBFC1C7F417C007B334C /* Porosity.hpp */,
BE60DBFD1C7F417C007B334C /* Porosity.hxx */,
BE51CDDD1CA4560100B66D2A /* TransientSource.cpp */,
BE51CDDE1CA4560100B66D2A /* TransientSource.hpp */,
BE51CDDF1CA4560100B66D2A /* TransientSource.hxx */,
);
path = GlobalVariationalOperatorInstances;
sourceTree = "<group>";
......@@ -7021,10 +7032,13 @@
BE669E1B1C7C9FDE00C521CA /* GradOnGradientBasedElasticityTensor.cpp */,
BE669E1C1C7C9FDE00C521CA /* GradOnGradientBasedElasticityTensor.hpp */,
BE669E1D1C7C9FDE00C521CA /* GradOnGradientBasedElasticityTensor.hxx */,
BE83AF2C1C733E15002C6FA3 /* HyperelasticLaw */,
BE51CDE21CA4560F00B66D2A /* TransientSource.cpp */,
BE51CDE31CA4560F00B66D2A /* TransientSource.hpp */,
BE51CDE41CA4560F00B66D2A /* TransientSource.hxx */,
BE60DBFF1C7F4222007B334C /* Porosity.cpp */,
BE60DC001C7F4222007B334C /* Porosity.hpp */,
BE60DC011C7F4222007B334C /* Porosity.hxx */,
BE83AF2C1C733E15002C6FA3 /* HyperelasticLaw */,
);
path = LocalVariationalOperatorInstances;
sourceTree = "<group>";
......@@ -12020,6 +12034,7 @@
BEF1832C1C6DE3B2008A6F1E /* Ale.cpp in Sources */,
BE83AF341C733EAD002C6FA3 /* StVenantKirchhoff.cpp in Sources */,
BE9BE1B51C96D02B0065BEFE /* SolidOnFluidMesh.cpp in Sources */,
BE51CDE51CA4560F00B66D2A /* TransientSource.cpp in Sources */,
BE5842131C6B8C920040C694 /* Mass.cpp in Sources */,
BEAEA8C71C71DFD30054FEF9 /* Porosity.cpp in Sources */,
BE669E1E1C7C9FDE00C521CA /* GradOnGradientBasedElasticityTensor.cpp in Sources */,
......@@ -12030,6 +12045,7 @@
BE669E1A1C7C9FCE00C521CA /* GradOnGradientBasedElasticityTensor.cpp in Sources */,
BEF183281C6DE3A3008A6F1E /* Ale.cpp in Sources */,
BE60DBFE1C7F417C007B334C /* Porosity.cpp in Sources */,
BE51CDE01CA4560100B66D2A /* TransientSource.cpp in Sources */,
BE83AF3F1C73669D002C6FA3 /* BulkSolid.cpp in Sources */,
);
runOnlyForDeploymentPostprocessing = 0;
......@@ -128,9 +128,9 @@ namespace HappyHeart
/*!
* \brief Solve the variational formulation for the current time step.
* \brief Set and solve the variational formulation for the current time step.
*
* This methid is intended to be called by the model's Forward() method.
* This method is intended to be called by the model's Forward() method.
*
* This covers the following steps:
* - Assembling matrices and vectors that ought to.
......
//! \file
//
//
// TransientSource.cpp
// HappyHeart
//
// Created by Sebastien Gilles on 03/02/15.
// Copyright (c) 2015 Inria. All rights reserved.
//
#include "ModelInstances/UnderDevelopment/Poromechanics/GlobalVariationalOperatorInstances/TransientSource.hpp"
namespace HappyHeart
{
namespace PoromechanicsNS
{
namespace GlobalVariationalOperatorNS
{
TransientSource::TransientSource(const FEltSpace& felt_space,
const Unknown& unknown,
const unsigned int geom_mesh_region_dimension,
const quadrature_rule_per_topology_type& quadrature_rule_per_topology,
const ParameterAtDof<ParameterNS::Type::scalar>& porosity,
DoComputeProcessorWiseLocal2Global do_consider_processor_wise_local_2_global,
const Parameter<ParameterNS::Type::vector>& source)
: parent(felt_space,
unknown,
geom_mesh_region_dimension,
quadrature_rule_per_topology,
AllocateGradientFEltPhi::no,
do_consider_processor_wise_local_2_global,
porosity,
source)
{ }
const std::string& TransientSource::ClassName()
{
static std::string name("Transient source with porosity");
return name;
}
} //namespace PoromechanicsNS
} // namespace GlobalVariationalOperatorNS
} // namespace HappyHeart
//! \file
//
//
// TransientSource.hpp
// HappyHeart
//
// Created by Sebastien Gilles on 03/02/15.
// Copyright (c) 2015 Inria. All rights reserved.
//
#ifndef HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_GLOBAL_VARIATIONAL_OPERATOR_INSTANCES_x_TRANSIENT_SOURCE_HPP_
# define HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_GLOBAL_VARIATIONAL_OPERATOR_INSTANCES_x_TRANSIENT_SOURCE_HPP_
# include "Parameters/Parameter.hpp"
# include "Parameters/ParameterAtDof.hpp"
# include "Operators/GlobalVariationalOperator/GlobalVariationalOperator.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/LocalVariationalOperatorInstances/TransientSource.hpp"
namespace HappyHeart
{
namespace PoromechanicsNS
{
namespace GlobalVariationalOperatorNS
{
/*!
* \brief Implementation of global TransientSource operator.
*
* \todo Improve the comment by writing its mathematical definition!
*/
class TransientSource final
: public GlobalVariationalOperator
<
TransientSource,
LocalVariationalOperatorNS::TransientSource
>
{
public:
//! Alias to unique pointer.
using const_unique_ptr = std::unique_ptr<const TransientSource>;
//! Returns the name of the operator.
static const std::string& ClassName();
//! Convenient alias to pinpoint the GlobalVariationalOperator parent.
using parent = GlobalVariationalOperator
<
TransientSource,
LocalVariationalOperatorNS::TransientSource
>;
//! Friendship to the parent class so that the CRTP can reach private methods defined below.
friend parent;
public:
/// \name Special members.
///@{
/*!
* \brief Constructor.
*
* \param[in] source Parameter describing the force applied.
* \param[in] felt_space Finite element space upon which the operator is defined.
* \param[in] unknown Unknown considered for this operator (might be scalar or vectorial).
* \param[in] geom_mesh_region_dimension Dimension of the geometric mesh region considered.
*/
explicit TransientSource(const FEltSpace& felt_space,
const Unknown& unknown,
unsigned int geom_mesh_region_dimension,
const quadrature_rule_per_topology_type& quadrature_rule_per_topology,
const ParameterAtDof<ParameterNS::Type::scalar>& porosity,
DoComputeProcessorWiseLocal2Global do_consider_processor_wise_local_2_global,
const Parameter<ParameterNS::Type::vector>& source);
//! Destructor.
~TransientSource() = default;
//! Copy constructor.
TransientSource(const TransientSource&) = default;
//! Move constructor.
TransientSource(TransientSource&&) = default;
//! Copy affectation.
TransientSource& operator=(const TransientSource&) = default;
//! Move affectation.
TransientSource& operator=(TransientSource&&) = default;
///@}
/*!
* \brief Assemble into one or several vectors.
*
* \tparam LinearAlgebraTupleT A tuple that may include \a GlobalVectorWithCoefficient objects.
*
* \param[in] global_vector_with_coeff_tuple List of global vectors into which the operator is
* assembled. These vectors are assumed to be already properly allocated.
* \param[in] domain Domain upon which the assembling takes place. Beware: if this domain is not a subset
* of the finite element space, assembling can only occur in a subset of the domain defined in the finite
* element space; if current \a domain is not a subset of finite element space one, assembling will occur
* upon the intersection of both.
* \param[in] time Time in seconds.
*/
template<class LinearAlgebraTupleT>
void Assemble(LinearAlgebraTupleT&& global_vector_with_coeff_tuple,
double time,
const Domain& domain = Domain()) const;
private:
/*!
* \brief Computes the additional arguments required by the AssembleWithVariadicArguments() method as
* a tuple.
*
* \param[in] local_operator Local operator in charge of the elementary computation. A mutable work
* variable is actually set in this call.
* \param[in] local_felt_space List of finite elements being considered; all those related to the
* same GeometricElt.
* \param[in] additional_arguments Additional arguments given to PerformElementaryCalculation().
* These arguments might need treatment before being given to ComputeEltArray: for instance if there
* is a GlobalVector that reflects a previous state, ComputeEltArray needs only the dofs that are
* relevant locally.
*/
void
SetComputeEltArrayArguments(const LocalFEltSpace& local_felt_space,
LocalVariationalOperator& local_operator,
std::tuple<double>&& additional_arguments) const;
};
} // namespace GlobalVariationalOperatorNS
} //namespace PoromechanicsNS
} // namespace HappyHeart
# include "ModelInstances/UnderDevelopment/Poromechanics/GlobalVariationalOperatorInstances/TransientSource.hxx"
#endif // HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_GLOBAL_VARIATIONAL_OPERATOR_INSTANCES_x_TRANSIENT_SOURCE_HPP_
//! \file
//
//
// TransientSource.hxx
// HappyHeart
//
// Created by Sebastien Gilles on 03/02/15.
// Copyright (c) 2015 Inria. All rights reserved.
//
#ifndef HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_GLOBAL_VARIATIONAL_OPERATOR_INSTANCES_x_TRANSIENT_SOURCE_HXX_
# define HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_GLOBAL_VARIATIONAL_OPERATOR_INSTANCES_x_TRANSIENT_SOURCE_HXX_
namespace HappyHeart
{
namespace PoromechanicsNS
{
namespace GlobalVariationalOperatorNS
{
template<class LinearAlgebraTupleT>
inline void TransientSource::Assemble(LinearAlgebraTupleT&& global_vector_with_coeff_tuple,
double time,
const Domain& domain) const
{
return parent::AssembleImpl(std::move(global_vector_with_coeff_tuple),
domain,
time);
}
inline void
TransientSource
::SetComputeEltArrayArguments(const LocalFEltSpace& local_felt_space,
LocalVariationalOperator& local_operator,
std::tuple<double>&& additional_arguments) const
{
static_cast<void>(local_felt_space);
local_operator.SetTime(std::get<0>(additional_arguments));
}
} // namespace GlobalVariationalOperatorNS
} //namespace PoromechanicsNS
} // namespace HappyHeart
#endif // HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_GLOBAL_VARIATIONAL_OPERATOR_INSTANCES_x_TRANSIENT_SOURCE_HXX_
......@@ -31,38 +31,86 @@ namespace HappyHeart
const Parameter<ParameterNS::Type::vector>& inlet_pressure,
DirichletBoundaryCondition::vector_shared_ptr&& boundary_condition_list)
: parent(mpi,
time_manager,
time_manager,
god_of_dof,
std::move(boundary_condition_list)),
porosity_parent(porosity),
inlet_pressure_(inlet_pressure)
inlet_pressure_(inlet_pressure),
fluid_velocity_numbering_subset_(god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)))
{ }
void VariationalFormulation::Perform()
{
decltype(auto) time_manager = parent::GetTimeManager();
const auto time = time_manager.GetTime();
decltype(auto) numbering_subset = GetFluidVelocityNumberingSubset();
auto& rhs = parent::GetNonCstSystemRhs(numbering_subset);
rhs.ZeroEntries(__FILE__, __LINE__);
GlobalVectorWithCoefficient vector(rhs, 1.);
GetInletPressureSourceOperator().Assemble(std::make_tuple(std::ref(vector)), time);
parent::ApplyEssentialBoundaryCondition<OperatorNS::Nature::nonlinear>(numbering_subset,
numbering_subset);
rhs.View(parent::MpiHappyHeart(), __FILE__, __LINE__);
// parent::SolveLinear<IsFactorized::no>(numbering_subset, numbering_subset);
// parent::WriteSolution(time_manager, numbering_subset);
}
void VariationalFormulation::AllocateMatricesAndVectors()
{
decltype(auto) numbering_subset = GetFluidVelocityNumberingSubset();
parent::AllocateSystemMatrix(numbering_subset, numbering_subset);
parent::AllocateSystemVector(numbering_subset);
}
void VariationalFormulation::SupplInit(const InputParameterList& input_parameter_data)
{
static_cast<void>(input_parameter_data);
{
// \todo #583: Temporary to make code work; quite stupid nonetheless!
const auto degree_of_exactness = DEFAULT_DEGREE_OF_EXACTNESS;
auto&& default_rule = DetermineDefaultQuadratureRule(degree_of_exactness);
this->GetNonCstQuadratureRulePerTopology().emplace_back(std::move(default_rule));
assert(this->GetNonCstQuadratureRulePerTopology().size() == 1ul);
}
decltype(auto) god_of_dof = parent::GetGodOfDof();
decltype(auto) unknown_manager = UnknownManager::GetInstance();
decltype(auto) fluid_velocity =
unknown_manager.GetUnknown(EnumUnderlyingType(UnknownIndex::fluid_velocity));
/*
const auto& inlet_border_felt_space = god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::inlet_border));
inlet_pressure_source_operator_ =
std::make_unique<GVO::TransientSource>(inlet_border_felt_space,
fluid_velocity,
mesh_dimension,
default_quadrature_rule_set,
DoComputeProcessorWiseLocal2Global::yes,
*inlet_pressure_);
const auto mesh_dimension = god_of_dof.GetGeometricMeshRegion().GetDimension();
const auto& default_quadrature_rule_set = this->GetQuadratureRulePerTopologyType();
const auto& inlet_border_felt_space =
god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::inlet_border));
namespace GVO = GlobalVariationalOperatorNS;
inlet_pressure_source_operator_ =
std::make_unique<GVO::TransientSource>(inlet_border_felt_space,
fluid_velocity,
mesh_dimension,
default_quadrature_rule_set,
GetPorosity(),
DoComputeProcessorWiseLocal2Global::yes,
GetInletPressure());
*/
}
......
......@@ -18,12 +18,13 @@
# include "Geometry/Domain.hpp"
# include "Operators/GlobalVariationalOperatorInstances/LinearForm/TransientSource.hpp"
# include "FormulationSolver/VariationalFormulation.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/InputParameterList.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/Crtp/Porosity.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/GlobalVariationalOperatorInstances/TransientSource.hpp"
namespace HappyHeart
......@@ -100,6 +101,24 @@ namespace HappyHeart
///@}
public:
/*!
* \todo #9 Put that Doxygen comment in common somewhere!
*
* \brief Set and solve the variational formulation for the current time step.
*
* This method is intended to be called by the model's Forward() method.
*
* This covers the following steps:
* - Assembling matrices and vectors that ought to.
* - Apply boundary conditions.
* - Call Petsc's solver to solve the system.
* - Write the solutions.
*/
void Perform();
private:
/// \name CRTP-required methods.
......@@ -135,12 +154,14 @@ namespace HappyHeart
///@{