Commit 22a3b67f authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#820 Poromechanics/implicit fluid step: introduce a differential enum class to...

#820 Poromechanics/implicit fluid step: introduce a differential enum class to reuse some of variational formulation in dH computation.
parent d30012d2
......@@ -664,6 +664,7 @@
BE5389FE1C897FE400D80749 /* Mpi.hxx in Headers */ = {isa = PBXBuildFile; fileRef = BE5389FC1C897FE400D80749 /* Mpi.hxx */; };
BE5657411D3CD7470091F063 /* Array.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BE56573F1D3CD7470091F063 /* Array.hpp */; };
BE5657421D3CD7470091F063 /* Array.hxx in Headers */ = {isa = PBXBuildFile; fileRef = BE5657401D3CD7470091F063 /* Array.hxx */; };
BE5657461D3CFF8F0091F063 /* Differential.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE5657431D3CFF8F0091F063 /* Differential.cpp */; };
BE5674161C722F3F00B2BBAD /* Policy.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BE5674121C722F3F00B2BBAD /* Policy.hpp */; };
BE5674171C722F3F00B2BBAD /* Policy.hxx in Headers */ = {isa = PBXBuildFile; fileRef = BE5674131C722F3F00B2BBAD /* Policy.hxx */; };
BE5674181C722F3F00B2BBAD /* VariationalFormulation.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BE5674141C722F3F00B2BBAD /* VariationalFormulation.hpp */; };
......@@ -4655,6 +4656,8 @@
BE5389FC1C897FE400D80749 /* Mpi.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = Mpi.hxx; sourceTree = "<group>"; };
BE56573F1D3CD7470091F063 /* Array.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = Array.hpp; sourceTree = "<group>"; };
BE5657401D3CD7470091F063 /* Array.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = Array.hxx; sourceTree = "<group>"; };
BE5657431D3CFF8F0091F063 /* Differential.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = Differential.cpp; sourceTree = "<group>"; };
BE5657441D3CFF8F0091F063 /* Differential.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = Differential.hpp; sourceTree = "<group>"; };
BE5674121C722F3F00B2BBAD /* Policy.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = Policy.hpp; sourceTree = "<group>"; };
BE5674131C722F3F00B2BBAD /* Policy.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = Policy.hxx; sourceTree = "<group>"; };
BE5674141C722F3F00B2BBAD /* VariationalFormulation.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = VariationalFormulation.hpp; sourceTree = "<group>"; };
......@@ -9367,6 +9370,8 @@
children = (
BEAE11F41D07072900967C19 /* VariationalFormulation.hpp */,
BEAE11F51D07072900967C19 /* VariationalFormulation.hxx */,
BE5657431D3CFF8F0091F063 /* Differential.cpp */,
BE5657441D3CFF8F0091F063 /* Differential.hpp */,
BEAE11EE1D07064A00967C19 /* InnerLoop */,
BEA78B951CE0925D00A6A185 /* GlobalVariationalOperatorInstances */,
BEA78B9C1CE0925D00A6A185 /* LocalVariationalOperatorInstances */,
......@@ -13350,6 +13355,7 @@
BEA78BB71CE092CA00A6A185 /* VariationalFormulation.cpp in Sources */,
BE82A89A1D0AA4EC00D3B579 /* Darcy.cpp in Sources */,
BE82A89C1D0AA4F200D3B579 /* Darcy.cpp in Sources */,
BE5657461D3CFF8F0091F063 /* Differential.cpp in Sources */,
BE60DBFE1C7F417C007B334C /* Porosity.cpp in Sources */,
BE9680051CB7F66E003F658B /* InternalFriction.cpp in Sources */,
BEAE0D251CB2AF1600CE333D /* ScalarDivVectorial.cpp in Sources */,
//! \file
//
//
// Differential.cpp
// HappyHeart
//
// Created by Sebastien Gilles on 18/07/16.
// Copyright © 2016 Inria. All rights reserved.
//
#include "Utilities/String/EmptyString.hpp"
#include "ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/Differential.hpp"
namespace HappyHeart
{
namespace PoromechanicsNS
{
namespace // anonymous
{
const std::string& Yes()
{
static std::string ret("differential_");
return ret;
}
} // namespace anonymous
const std::string& DifferentialPreffix(differential is_differential)
{
switch (is_differential)
{
case differential::yes:
return Yes();
case differential::no:
return Utilities::EmptyString();
}
}
} // namespace PoromechanicsNS
} // namespace HappyHeart
//! \file
//
//
// Differential.hpp
// HappyHeart
//
// Created by Sebastien Gilles on 18/07/16.
// Copyright © 2016 Inria. All rights reserved.
//
#ifndef HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_IMPLICIT_STEP_FLUID_x_DIFFERENTIAL_HPP_
# define HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_IMPLICIT_STEP_FLUID_x_DIFFERENTIAL_HPP_
# include <memory>
# include <vector>
namespace HappyHeart
{
namespace PoromechanicsNS
{
//! Enum class to choose whether it is run in differential mode or not.
enum class differential { no, yes };
//! Returns preffix "differential_" if differential is yes.
const std::string& DifferentialPreffix(differential is_differential);
} // namespace PoromechanicsNS
} // namespace HappyHeart
#endif // HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_IMPLICIT_STEP_FLUID_x_DIFFERENTIAL_HPP_
......@@ -33,6 +33,7 @@
# include "ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/GlobalVariationalOperatorInstances/T21.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/GlobalVariationalOperatorInstances/Darcy.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/GlobalVariationalOperatorInstances/HybridVector.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/Differential.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/Crtp/VariableHolder.hpp"
......@@ -44,14 +45,14 @@ namespace HappyHeart
{
namespace ImplicitStepFluidNS
{
enum class DoApplyBoundaryCondition { yes, no };
enum class DoRecomputeBlock { yes, no };
namespace InnerLoopNS
{
......@@ -173,7 +174,7 @@ namespace HappyHeart
*
* \return Norm obtained.
*/
double Perform(unsigned int iteration);
double Perform(differential is_differential, unsigned int iteration);
//! Compute the matrix with an updated T21 block (to use after the Newton converged).
const GlobalMatrix& ComputeMatrixWithUpdatedT21() noexcept;
......@@ -546,7 +547,7 @@ namespace HappyHeart
//! Compute rhs.
void ComputeRhs();
void ComputeRhs(differential is_differential);
//! Compute T11.
void ComputeT11(Wrappers::Petsc::DoReuseMatrix do_reuse_matrix,
......@@ -836,7 +837,7 @@ namespace HappyHeart
} // namespace ImplicitStepFluidNS
} // namespace PoromechanicsNS
......
......@@ -29,7 +29,9 @@ namespace HappyHeart
template<class HyperelasticLawT>
double VariationalFormulation<HyperelasticLawT>::Perform(unsigned int iteration)
double VariationalFormulation<HyperelasticLawT>
::Perform(differential is_differential,
unsigned int iteration)
{
SetIndexInnerLoop(iteration);
decltype(auto) time_manager = parent::GetTimeManager();
......@@ -71,21 +73,30 @@ namespace HappyHeart
system_matrix.ZeroEntries(__FILE__, __LINE__);
ComputeMatrix(do_reuse_matrix, BlocksToRecompute::all);
ComputeRhs(); // must occur after T11, as it needs the t11 matrix.
ComputeRhs(is_differential); // must occur after T11, as it needs the t11 matrix.
system_matrix.Assembly(__FILE__, __LINE__);
system_matrix.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset) + "/full_matrix_" + GetIterationTag() + ".hhdata",
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "full_matrix_" + GetIterationTag() + ".hhdata",
__FILE__, __LINE__);
system_matrix.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset) + "/full_matrix_" + GetIterationTag() + ".m",
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "full_matrix_" + GetIterationTag() + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
system_matrix.ViewBinary(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset) + "/full_matrix_" + GetIterationTag() + ".binary.hhdata",
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "full_matrix_" + GetIterationTag() + ".binary.hhdata",
__FILE__, __LINE__);
......@@ -96,7 +107,10 @@ namespace HappyHeart
// \todo #820 TMP For dev: write the solution in one file, aggregating from all ranks.
auto& solution = parent::GetNonCstSystemSolution(numbering_subset);
solution.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset) + "/correction_" + GetIterationTag() + ".m",
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "correction_" + GetIterationTag() + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
......@@ -111,7 +125,10 @@ namespace HappyHeart
corrected_values.UpdateGhosts(__FILE__, __LINE__);
corrected_values.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset) + "/mixt_variables_" + GetIterationTag() + ".m",
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "mixt_variables_" + GetIterationTag() + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
......@@ -243,7 +260,7 @@ namespace HappyHeart
template<class HyperelasticLawT>
void VariationalFormulation<HyperelasticLawT>::ComputeRhs()
void VariationalFormulation<HyperelasticLawT>::ComputeRhs(differential is_differential)
{
decltype(auto) god_of_dof = parent::GetGodOfDof();
decltype(auto) bc_manager = DirichletBoundaryConditionManager::GetInstance();
......@@ -258,7 +275,10 @@ namespace HappyHeart
// \todo #820 I assume here Darcy is the first thing assembled into rhs. Anyway this line will
// disappear once model is validated.
rhs.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset) + "/darcy_" + GetIterationTag() + ".m",
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "darcy_" + GetIterationTag() + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
......@@ -268,7 +288,10 @@ namespace HappyHeart
auto& contribution = GetNonCstWorkRhsLike(); // \todo #820 Intermediate step for dev purpose only.
contribution.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset) + "/pressure_contrib_to_rhs_" + GetIterationTag() + ".m",
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "pressure_contrib_to_rhs_" + GetIterationTag() + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
......@@ -319,7 +342,10 @@ namespace HappyHeart
contribution.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset) + "/fluid_mass_contrib_to_rhs_" + GetIterationTag() + ".m",
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "fluid_mass_contrib_to_rhs_" + GetIterationTag() + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
......@@ -339,7 +365,10 @@ namespace HappyHeart
rhs.UpdateGhosts(__FILE__, __LINE__);
rhs.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset) + "/rhs_" + GetIterationTag() + ".m",
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "rhs_" + GetIterationTag() + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
......@@ -400,7 +429,6 @@ namespace HappyHeart
return parent::GetSystemMatrix(numbering_subset, numbering_subset);
}
} // namespace InnerLoopNS
......
......@@ -35,9 +35,8 @@ namespace HappyHeart
namespace ImplicitStepFluidNS
{
template<class HyperelasticLawT>
class VariationalFormulation
: public HappyHeart::VariationalFormulation
......@@ -130,8 +129,10 @@ namespace HappyHeart
* - Apply boundary conditions.
* - Call Petsc's solver to solve the system.
* - Write the solutions.
*
* \param[in] is_differential If yes, run the differential of the variational formulation.
*/
void Perform();
void Perform(differential is_differential);
//! Constant accessor to the darcy vector after inner loop.
const GlobalVector& GetDarcyVector() const noexcept;
......
......@@ -58,7 +58,7 @@ namespace HappyHeart
template<class HyperelasticLawT>
void VariationalFormulation<HyperelasticLawT>::Perform()
void VariationalFormulation<HyperelasticLawT>::Perform(differential is_differential)
{
double current_norm = 1.e20;
......@@ -72,7 +72,7 @@ namespace HappyHeart
// Newton loop (\todo #820 To replace with Petsc Newton...)
do
{
current_norm = inner_varf.Perform(counter++);
current_norm = inner_varf.Perform(is_differential, counter++);
std::cout << "Norm after inner loop = " << current_norm << std::endl;
} while (current_norm > 1.e-12 && counter < 10u); // \todo #820 Replace by proper stop condition from input file!
......@@ -90,7 +90,11 @@ namespace HappyHeart
GetNonCstMonolithicMatrixAfterInnerLoop().Copy(updated_matrix, __FILE__, __LINE__);
updated_matrix.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset) + "/matrix_after_FP_mixt_with_dirichlet_time_" + std::to_string(parent::GetTimeManager().NtimeModified()) + ".m",
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "matrix_after_FP_mixt_with_dirichlet_time_"
+ std::to_string(parent::GetTimeManager().NtimeModified()) + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
......@@ -103,7 +107,11 @@ namespace HappyHeart
// \todo #820 Not exactly equal to Freefem; investigate if one term is missing!
updated_matrix.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset) + "/matrix_after_FP_mixt_without_dirichlet_time_" + std::to_string(parent::GetTimeManager().NtimeModified()) + ".m",
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "matrix_after_FP_mixt_without_dirichlet_time_"
+ std::to_string(parent::GetTimeManager().NtimeModified()) + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
......@@ -114,7 +122,11 @@ namespace HappyHeart
inner_varf.ComputeDarcy(-1., IsFullDarcy::no, darcy_vector);
darcy_vector.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset) + "/darcy_vector_after_FP_mixt_time_" + std::to_string(parent::GetTimeManager().NtimeModified()) + ".m",
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "darcy_vector_after_FP_mixt_time_"
+ std::to_string(parent::GetTimeManager().NtimeModified()) + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
......
......@@ -48,7 +48,7 @@ namespace HappyHeart
}
auto& implicit_step_fluid_varf = GetNonCstImplicitStepFluidVarf();
implicit_step_fluid_varf.Perform();
implicit_step_fluid_varf.Perform(differential::no);
UpdateSolidVector();
......
......@@ -16,6 +16,7 @@ poromechanics_src = Split('''
./GlobalVariationalOperatorInstances/Porosity.cpp
./GlobalVariationalOperatorInstances/ScalarDivVectorial.cpp
./HyperelasticLaw/StVenantKirchhoff.cpp
./ImplicitStepFluid/Differential.cpp
./ImplicitStepFluid/GlobalVariationalOperatorInstances/Darcy.cpp
./ImplicitStepFluid/GlobalVariationalOperatorInstances/HybridVector.cpp
./ImplicitStepFluid/GlobalVariationalOperatorInstances/T33.cpp
......@@ -37,7 +38,7 @@ poromechanics_src = Split('''
./Porosity/VariationalFormulation.cpp
./Private/SolidOnFluidMesh.cpp
./VariableHolder.cpp
''')
''')
poromechanics_lib = custom.HappyHeartLibrary(env, 'happy_heart_poromechanics', poromechanics_src, [suppl_link_libraries])
......
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