Commit 05487282 authored by thomaskummer's avatar thomaskummer
Browse files

restart

parent 913ef769
......@@ -528,6 +528,7 @@ void ExporterHDF5<MeshType>::import (const Real& Tstart, const Real& dt)
template<typename MeshType>
void ExporterHDF5<MeshType>::import (const Real& time)
{
std::cout << time << std::endl;
if ( M_HDF5.get() == 0)
{
M_HDF5.reset (new hdf5_Type (this->M_dataVector.begin()->storedArrayPtr()->comm() ) );
......
......@@ -29,7 +29,7 @@ hdf5_sheets_name = sheets
[./boundary_conditions]
list = 'EndocardiumLV Epicardium EndocardiumRV Septum Support Base SeptumEdge'
list = 'EndocardiumLV Epicardium EndocardiumRV Septum Support Base'
[./EndocardiumLV]
type = Essential
......
......@@ -141,7 +141,8 @@ template<typename Mesh> inline void importVectorField ( boost::shared_ptr<Vecto
const std::string& fieldName,
boost::shared_ptr< Mesh > localMesh ,
const std::string& postDir = "./",
const std::string& polynomialDegree = "P1" )
const std::string& polynomialDegree = "P1",
const std::string& timestep = "00000" )
{
typedef Mesh mesh_Type;
typedef ExporterData<mesh_Type> exporterData_Type;
......@@ -153,7 +154,7 @@ template<typename Mesh> inline void importVectorField ( boost::shared_ptr<Vecto
boost::shared_ptr<Epetra_Comm> comm ( new Epetra_MpiComm (MPI_COMM_WORLD) );
boost::shared_ptr<FESpace< mesh_Type, MapEpetra > > feSpace ( new FESpace< mesh_Type, MapEpetra > ( localMesh, polynomialDegree, 3, comm ) );
exporterData_Type impData (exporterData_Type::VectorField, fieldName + ".00000", feSpace,
exporterData_Type impData (exporterData_Type::VectorField, fieldName + "." + timestep, feSpace,
vector, UInt (0), exporterData_Type::UnsteadyRegime);
filterPtr_Type importer ( new hdf5Filter_Type() );
......
......@@ -18,7 +18,7 @@ ADD_SUBDIRECTORIES(
# linearizedGHOrotation
# linearizedGHOventricle
# prestress
# example_EMMaterial
example_EMMaterial
example_HeartOharaRudy
# example_passiveEMMaterial
example_EMCoupling
......
......@@ -129,20 +129,11 @@ monodomain_xml_file = ParamList.xml
[../boundary_conditions]
#list = 'Support Base'
list = 'SupportNormal SupportRobin BaseNormal BaseRobin'
listVariableBC = 'LVEndo RVEndo RVSeptum RVRadii'
numPreloadSteps = 1
numPreloadSteps = 10
[./Support]
type = Essential
flag = 49
mode = Full
component = 3
function = '[0.0, 0.0, 0.0]'
[../SupportNormal]
[./SupportNormal]
type = Essential
flag = 49
mode = Component
......@@ -164,13 +155,6 @@ monodomain_xml_file = ParamList.xml
[../]
[../Base]
type = Essential
flag = 48
mode = Full
component = 3
function = '[0.0, 0.0, 0.0]'
[../BaseNormal]
type = Essential
flag = 48
......@@ -215,7 +199,7 @@ monodomain_xml_file = ParamList.xml
pPerturbationFe = 1e-2
pPerturbationCirc = 1e-3
dpMax = 5.0
couplingError = 1e-8
couplingError = 1e-7
couplingJFeIter = 20
couplingJFeSubStart = 2
couplingJFeSubIter = 5
......
......@@ -314,8 +314,7 @@ int main (int argc, char** argv)
ExporterHDF5< RegionMesh <LinearTetra> > activationTimeExporter;
activationTimeExporter.setMeshProcId (solver.localMeshPtr(), solver.comm()->MyPID() );
activationTimeExporter.addVariable (ExporterData<mesh_Type>::ScalarField, "Activation Time",
solver.electroSolverPtr()->feSpacePtr(), activationTimeVector, UInt (0) );
activationTimeExporter.addVariable (ExporterData<mesh_Type>::ScalarField, "Activation Time", solver.electroSolverPtr()->feSpacePtr(), activationTimeVector, UInt (0) );
activationTimeExporter.setPrefix ("ActivationTime");
activationTimeExporter.setPostDir (problemFolder);
......@@ -332,10 +331,9 @@ int main (int argc, char** argv)
//============================================//
const std::string circulationInputFile = command_line.follow ("inputfile", 2, "-cif", "--cifile");
const std::string circulationOutputFile = command_line.follow ("solution.txt", 2, "-cof", "--cofile");
const std::string circulationOutputFile = command_line.follow ( (problemFolder + "solution.txt").c_str(), 2, "-cof", "--cofile");
Circulation circulationSolver( circulationInputFile );
if ( 0 == comm->MyPID() ) circulationSolver.exportSolution( circulationOutputFile );
// Flow rate between two vertices
auto Q = [&circulationSolver] (const std::string& N1, const std::string& N2) { return circulationSolver.solution ( std::vector<std::string> {N1, N2} ); };
......@@ -455,7 +453,8 @@ int main (int argc, char** argv)
UInt iter;
Real t (0);
auto printCoupling = [&] ( std::string label ) { if ( 0 == comm->MyPID() ) {
auto printCoupling = [&] ( std::string label ) { if ( 0 == comm->MyPID() )
{
std::cout << "\n****************************** Coupling: " << label << " *******************************";
std::cout << "\nNewton iteration nr. " << iter << " at time " << t;
std::cout << "\nLV - Pressure: " << bcValues[0];
......@@ -472,35 +471,90 @@ int main (int argc, char** argv)
std::cout << "\n****************************** Coupling: " << label << " *******************************\n\n"; }
};
auto pipeToString = [] ( const char* command )
{
FILE* file = popen( command, "r" ) ;
std::ostringstream stm ;
char line[6] ;
fgets( line, 6, file );
stm << line;
pclose(file) ;
return stm.str() ;
};
//============================================//
// Preload
// Load restart file
//============================================//
const int preloadSteps = dataFile ( "solid/boundary_conditions/numPreloadSteps", 0);
auto preloadPressure = [] (std::vector<double> p, const int& step, const int& steps)
std::string restartInput = command_line.follow ("noRestart", 2, "-r", "--restart");
const bool restart ( restartInput != "noRestart" );
if ( restart )
{
for (auto& i : p) {i *= double(step) / double(steps);}
return p;
};
// Get most recent restart index
if ( restartInput == "." ) restartInput = pipeToString( ("grep Iteration " + problemFolder + "MechanicalSolution.xmf | tail -n 1 | awk '{print $5}'").c_str() );
// Set time variable
const unsigned int nIter = (std::stoi(restartInput) - 1) / saveIter;
t = nIter * dt_mechanics;
std::cout << nIter << " " << restartInput << " " << std::stoi(restartInput) << " " << saveIter << std::endl;
activationTimeExporter.setTimeIndex(nIter * saveIter + 1);
solver.setTimeIndex(nIter * saveIter + 1);
// FE
const std::string restartDir = command_line.follow (problemFolder.c_str(), 2, "-rd", "--restartDir");
ElectrophysiologyUtility::importVectorField ( solver.structuralOperatorPtr() -> displacementPtr(), "MechanicalSolution" , "displacement", solver.localMeshPtr(), restartDir, "P2", restartInput );
//solver.electroSolverPtr() -> importSolution ("ElectroSolution", problemFolder, t);
ElectrophysiologyUtility::importVectorField (solver.electroSolverPtr()->potentialPtr(), "ElectroSolution" , "Variable0", solver.localMeshPtr(), restartDir, "P2", restartInput );
//ElectrophysiologyUtility::importVectorField (solver.activationModelPtr() -> fiberActivationPtr(), "ActivationSolution" , "Activation", solver.localMeshPtr(), restartDir, "P2", restartInput );
// Circulation
circulationSolver.restartFromFile ( circulationOutputFile , nIter );
// Coupling boundary conditions
bcValues = { p ( "lv" ) , p ( "rv") };
}
solver.saveSolution (-1.0);
//============================================//
// Preload
//============================================//
for (int i (1); i <= preloadSteps; i++)
if ( ! restart )
{
if ( 0 == comm->MyPID() )
const int preloadSteps = dataFile ( "solid/boundary_conditions/numPreloadSteps", 0);
auto preloadPressure = [] (std::vector<double> p, const int& step, const int& steps)
{
std::cout << "\n*********************";
std::cout << "\nPreload step: " << i << " / " << preloadSteps;
std::cout << "\n*********************\n";
}
for (auto& i : p) {i *= double(step) / double(steps);}
return p;
};
solver.saveSolution (-1.0);
activationTimeExporter.postProcess(-1.0);
// Update pressure b.c.
modifyFeBC(preloadPressure(bcValues, i, preloadSteps));
for (int i (1); i <= preloadSteps; i++)
{
if ( 0 == comm->MyPID() )
{
std::cout << "\n*********************";
std::cout << "\nPreload step: " << i << " / " << preloadSteps;
std::cout << "\n*********************\n";
}
// Update pressure b.c.
modifyFeBC(preloadPressure(bcValues, i, preloadSteps));
// Solve mechanics
solver.bcInterfacePtr() -> updatePhysicalSolverVariables();
solver.solveMechanics();
// Solve mechanics
solver.bcInterfacePtr() -> updatePhysicalSolverVariables();
solver.solveMechanics();
}
}
......@@ -512,6 +566,8 @@ int main (int argc, char** argv)
VFe[1] = RV.volume(disp, dETFESpace, 1);
VCirc = VFe;
printCoupling("Initial values");
auto perturbedPressure = [] (std::vector<double> p, const double& dp)
{
for (auto& i : p) {i += dp;}
......@@ -525,7 +581,9 @@ int main (int argc, char** argv)
};
solver.saveSolution (t);
activationTimeExporter.postProcess(t);
if ( 0 == comm->MyPID() ) circulationSolver.exportSolution( circulationOutputFile , restart );
for (int k (1); k <= maxiter; k++)
{
if ( 0 == comm->MyPID() )
......
......@@ -10,10 +10,10 @@
<Parameter name="elementsOrder" type="string" value="P1" />
<Parameter name="ionic_model" type="string" value="AlievPanfilov" />
<Parameter name="solid_mesh_name" type="string" value="cube4.mesh" />
<Parameter name="solid_mesh_path" type="string" value="/Users/srossi/LifeV/meshes/" />
<Parameter name="solid_mesh_path" type="string" value="../../data/mesh/" />
<Parameter name="solid_fiber_file" type="string" value="idealHeart_structure_fiber" />
<Parameter name="mesh_name" type="string" value="cube4.mesh" />
<Parameter name="mesh_path" type="string" value="/Users/srossi/LifeV/meshes/" />
<Parameter name="mesh_path" type="string" value="/Users/thomas/LIFE5/lifev-em/lifev/em/data/mesh/" />
<Parameter name="fiber_file" type="string" value="idealHeart_electro_fiber" />
<Parameter name="fiber_X" type="double" value="0.707106781187" />
<Parameter name="fiber_Y" type="double" value="0.707106781187" />
......
......@@ -55,7 +55,7 @@ save = 1
# InverseViscosity
# ActiveForceCoefficient
# ChemicalThreshold
ActivationModel = ActiveStrainRossi14
ActivationModel = ActiveStressRossi14
CalciumIndex = 0
ActiveForceCoefficient = -4.0
......@@ -102,8 +102,8 @@ lawType = nonlinear
#ANH = Active (Strain) Neo Hookean
#SimpleActiveStress = Simple Fiber Active Stress
EMPassiveMaterialType = PV
EMActiveStressMaterialType = ANH #SFAS
EMPassiveMaterialType = PHO
EMActiveStressMaterialType = SimpleActiveStress #SFAS
Tmax = 497000
mu = 4960
......@@ -182,8 +182,8 @@ BDF_order = 2
[../space_discretization]
mesh_type = .mesh
mesh_dir = /usr/scratch/srossi/meshes/
mesh_file = CoarseSymmCube.mesh
mesh_dir = ../../data/mesh/
mesh_file = Cube.mesh
order = P1
......
......@@ -326,6 +326,9 @@ public:
}
void saveSolution (Real time);
void setTimeIndex (const UInt& time);
void closeExporters();
......@@ -632,6 +635,15 @@ EMSolver<Mesh, ElectroSolver>::setupExporters (std::string problemFolder,
UInt (0) );
}
template<typename Mesh , typename ElectroSolver>
void
EMSolver<Mesh, ElectroSolver>::setTimeIndex (const UInt& time)
{
M_electroExporterPtr -> setTimeIndex (time);
M_activationExporterPtr -> setTimeIndex (time);
M_mechanicsExporterPtr -> setTimeIndex (time);
}
template<typename Mesh , typename ElectroSolver>
void
EMSolver<Mesh, ElectroSolver>::saveSolution (Real time)
......
......@@ -149,7 +149,12 @@ public:
void updateTimeVar(const double& dt)
{
M_time += dt;
setTime(M_time + dt);
}
void setTime(const double& time)
{
M_time = time;
}
void updatePrevSolVectors()
......@@ -198,17 +203,46 @@ public:
std::cout << "t = " << M_time << ":\t\t" << "rel. error = " << relative_error << std::endl;
}
void exportSolution(const std::string& filename) const
void exportSolution(const std::string& filename, const bool& restart = false) const
{
VectorStdDouble exportRow ({M_time});
VectorStdDouble u ( eigenToStd(M_u) );
exportRow.insert(exportRow.end(), u.begin(), u.end());
CirculationIO exporter;
const bool append = (M_time == 0 ? false : true);
const bool append = (M_time == 0 && restart == false ? false : true);
exporter.exportVector(filename, exportRow, append);
}
void restartFromFile(const std::string& filename, const unsigned int& timestep)
{
CirculationIO importer;
std::vector<double> u ( importer.importVector(filename, timestep) );
std::vector<double> uPrev0 ( importer.importVector(filename, (timestep < 1 ? 0 : (timestep - 1))) );
std::vector<double> uPrev1 ( importer.importVector(filename, (timestep < 2 ? 0 : (timestep - 2))) );
M_u = stdToEigen(u).tail( M_u.size() );
M_uPrev0 = ( stdToEigen(uPrev0) ).tail( M_u.size() );
M_uPrev1 = ( stdToEigen(uPrev1) ).tail( M_u.size() );
M_time = u[0];
DofHandler dofh (M_gv);
for ( auto& element : M_gv.elements() )
{
const unsigned int elementIdx = dofh( element );
const unsigned int vertex1Idx = dofh( element->node(0) );
const unsigned int vertex2Idx = dofh( element->node(1) );
std::vector<double> u {M_u(elementIdx) , 0.0, 0.0 };
if ( vertex1Idx < dofh.size() ) u[1] = M_u(vertex1Idx);
if ( vertex2Idx < dofh.size() ) u[2] = M_u(vertex2Idx);
element -> initRestart ( u , M_time );
}
}
private:
......
......@@ -32,6 +32,7 @@ public:
const std::string node(const unsigned int idx) const {return M_nodes[idx];} // todo: return pointer to Vertex object
const std::vector<std::string> nodes() const {return M_nodes;}
const double init() const {return M_init;}
virtual void initRestart(const std::vector<double>& u, const double& time){}
virtual const double rhs(const std::vector<double>& u, const double& time) = 0;
virtual const double lhs(const unsigned int& idx, const std::vector<double>& u, const double& time) = 0; // u containing (Q p1 p2 R)
......@@ -311,6 +312,14 @@ public:
}
}
virtual void initRestart(const std::vector<double>& u, const double& time)
{
M_time = time;
M_D = M_param[0] * u[0] / ( u[1] - u[2] );
M_Dprev0 = M_D;
M_Dprev1 = M_D;
}
protected:
......@@ -353,6 +362,15 @@ public:
break;
};
}
virtual void initRestart(const std::vector<double>& u, const double& time)
{
M_time = time;
M_D = M_param[0] * u[0] * u[1] / ( u[1] - u[2] );
M_Dprev0 = M_D;
M_Dprev1 = M_D;
}
};
......
......@@ -183,6 +183,40 @@ public:
else std::cout << "Unable to open file";
}
const std::vector<double> importVector(const std::string& filename, const unsigned int& row) const
{
std::vector<double> vec;
std::ifstream file (filename.c_str());
if (file.is_open())
{
std::string line;
unsigned int i = 0;
while ( getline (file, line) )
{
std::istringstream iss(line);
std::string column;
if ( i == row )
{
while(iss >> column)
{
std::istringstream entryStr(column);
double entry;
entryStr >> entry;
vec.push_back ( entry );
}
}
++i;
}
file.close();
return vec;
}
else std::cout << "Unable to open file";
}
private:
......
......@@ -45,10 +45,8 @@ PassiveHolzapfelOgden<Mesh>::PassiveHolzapfelOgden() :
this -> M_materialFunctionList[3].reset ( new MaterialFunctions::dIsotropicExponential<Mesh> (3330.0, 9.242) );
this -> M_materialFunctionList[4].reset ( new MaterialFunctions::AnisotropicExponential<Mesh> (185350, 15.972) );
this -> M_materialFunctionList[5].reset ( new MaterialFunctions::dAnisotropicExponential<Mesh> (185350, 15.972) );
this -> M_materialFunctionList[6].reset (
new MaterialFunctions::AnisotropicExponential<Mesh> ( 25640.0, 10.446, MaterialFunctions::AnisotropicExponential<Mesh>::Sheets ) );
this -> M_materialFunctionList[7].reset (
new MaterialFunctions::dAnisotropicExponential<Mesh> ( 25640.0, 10.446, MaterialFunctions::AnisotropicExponential<Mesh>::Sheets ) );
this -> M_materialFunctionList[6].reset ( new MaterialFunctions::AnisotropicExponential<Mesh> ( 25640.0, 10.446, MaterialFunctions::AnisotropicExponential<Mesh>::Sheets ) );
this -> M_materialFunctionList[7].reset ( new MaterialFunctions::dAnisotropicExponential<Mesh> ( 25640.0, 10.446, MaterialFunctions::AnisotropicExponential<Mesh>::Sheets ) );
this -> M_materialFunctionList[8].reset ( new MaterialFunctions::ShearExponential<Mesh> (4170.0, 11.602) );
this -> M_materialFunctionList[9].reset ( new MaterialFunctions::dShearExponential<Mesh> (4170.0, 11.602) );
......
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