Commit ade5032a authored by thomaskummer's avatar thomaskummer
Browse files

intrepid lifev update

parent 453ac5f7
......@@ -457,16 +457,6 @@ int main (int argc, char** argv)
std::cout << "\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() ;
};
//============================================//
......@@ -478,64 +468,13 @@ int main (int argc, char** argv)
if ( restart )
{
const std::string restartDir = command_line.follow (problemFolder.c_str(), 2, "-rd", "--restartDir");
Real dtExport = 10.;
// Get most recent restart index
if ( restartInput == "." )
{
restartInput = pipeToString( ("tail -n 1 " + restartDir + "solution.dat | awk -F '[. ]' '{print $1 \".\" $2}' | awk '{printf \"%05g\", int($1*1000/" + std::to_string(dtExport) + ") + 1}'").c_str() );
}
// Set time variable
const unsigned int restartInputStr = std::stoi(restartInput);
const unsigned int nIter = (restartInputStr - 1) * dtExport / heartSolver.data().dt_mechanics();
t = nIter * heartSolver.data().dt_mechanics();
if ( 0 == comm->MyPID() )
{
std::cout << "\nLoad from restart: " << restartInput << ", nIterCirculation = " << nIter << ", time = " << t << std::endl;
}
heartSolver.restart(command_line);
// Set time exporter time index
solver.setTimeIndex(restartInputStr + 1);
//solver.importHdf5();
if ( 0 == comm->MyPID() )
{
std::cout << "\nLoad from restart: " << restartInput << ", nIterCirculation = " << nIter << ", time = " << t << std::endl;
}
// Load restart solutions from output files
std::string polynomialDegree = dataFile ( "solid/space_discretization/order", "P1");
if ( 0 == comm->MyPID() )
{
std::cout << "\nLoad from restart: " << restartInput << ", nIterCirculation = " << nIter << ", time = " << t << std::endl;
}
ElectrophysiologyUtility::importVectorField ( solver.structuralOperatorPtr() -> displacementPtr(), "MechanicalSolution" , "displacement", solver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
if ( 0 == comm->MyPID() )
{
std::cout << "\nLoad from restart: " << restartInput << ", nIterCirculation = " << nIter << ", time = " << t << std::endl;
}
for ( unsigned int i = 0; i < solver.electroSolverPtr()->globalSolution().size() ; ++i )
{
ElectrophysiologyUtility::importScalarField (solver.electroSolverPtr()->globalSolution().at(i), "ElectroSolution" , ("Variable" + std::to_string(i)), solver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
}
if ( 0 == comm->MyPID() )
{
std::cout << "\nLoad from restart: " << restartInput << ", nIterCirculation = " << nIter << ", time = " << t << std::endl;
}
ElectrophysiologyUtility::importScalarField (solver.activationModelPtr() -> fiberActivationPtr(), "ActivationSolution" , "Activation", solver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
//ElectrophysiologyUtility::importScalarField (solver.activationTimePtr(), "ActivationTimeSolution" , "Activation Time", solver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
if ( 0 == comm->MyPID() )
{
std::cout << "\nLoad from restart: " << restartInput << ", nIterCirculation = " << nIter << ", time = " << t << std::endl;
}
circulationSolver.restartFromFile ( restartDir + "solution.dat" , nIter );
// Set boundary mechanics conditions
bcValues = { p ( "lv" ) , p ( "rv" ) };
bcValuesPre = { p ( "lv" ) , p ( "rv" ) };
modifyFeBC(bcValues);
}
......@@ -546,55 +485,6 @@ int main (int argc, char** argv)
if ( ! restart )
{
heartSolver.preload(modifyFeBC, bcValues);
solver.structuralOperatorPtr() -> data() -> dataTime() -> setTime(0.0);
auto preloadPressure = [] (std::vector<double> p, const int& step, const int& steps)
{
for (auto& i : p) {i *= double(step) / double(steps);}
return p;
};
LifeChrono chronoSave;
chronoSave.start();
solver.saveSolution (-1.0);
if ( 0 == comm->MyPID() )
{
std::cout << "\n*****************************************************************";
std::cout << "\nData stored in " << chronoSave.diff() << " s";
std::cout << "\n*****************************************************************\n";
}
LifeChrono chronoPreload;
chronoPreload.start();
for (int i (1); i <= heartSolver.data().preloadSteps(); i++)
{
if ( 0 == comm->MyPID() )
{
std::cout << "\n*****************************************************************";
std::cout << "\nPreload step: " << i << " / " << heartSolver.data().preloadSteps();
std::cout << "\n*****************************************************************\n";
}
// Update pressure b.c.
modifyFeBC(preloadPressure(bcValues, i, heartSolver.data().preloadSteps() ));
// Solve mechanics
solver.bcInterfacePtr() -> updatePhysicalSolverVariables();
solver.solveMechanics();
//solver.saveSolution (i-1);
}
if ( 0 == comm->MyPID() )
{
std::cout << "\n*****************************************************************";
std::cout << "\nPreload done in: " << chronoPreload.diff();
std::cout << "\n*****************************************************************\n";
}
}
......
......@@ -215,6 +215,64 @@ public:
}
void restart(const GetPot& command_line)
{
const std::string restartDir = command_line.follow (problemFolder.c_str(), 2, "-rd", "--restartDir");
Real dtExport = 10.;
// Get most recent restart index
if ( restartInput == "." )
{
restartInput = pipeToString( ("tail -n 1 " + restartDir + "solution.dat | awk -F '[. ]' '{print $1 \".\" $2}' | awk '{printf \"%05g\", int($1*1000/" + std::to_string(dtExport) + ") + 1}'").c_str() );
}
// Set time variable
const unsigned int restartInputStr = std::stoi(restartInput);
const unsigned int nIter = (restartInputStr - 1) * dtExport / data().dt_mechanics();
t = nIter * data().dt_mechanics();
if ( 0 == M_emSolver.comm()->MyPID() )
{
std::cout << "\nLoad from restart: " << restartInput << ", nIterCirculation = " << nIter << ", time = " << t << std::endl;
}
// Set time exporter time index
M_emSolver.setTimeIndex(restartInputStr + 1);
//solver.importHdf5();
if ( 0 == M_emSolver.comm()->MyPID() )
{
std::cout << "\nLoad from restart: " << restartInput << ", nIterCirculation = " << nIter << ", time = " << t << std::endl;
}
// Load restart solutions from output files
std::string polynomialDegree = data().datafile() ( "solid/space_discretization/order", "P1");
if ( 0 == M_emSolver.comm()->MyPID() )
{
std::cout << "\nLoad from restart: " << restartInput << ", nIterCirculation = " << nIter << ", time = " << t << std::endl;
}
ElectrophysiologyUtility::importVectorField ( M_emSolver.structuralOperatorPtr() -> displacementPtr(), "MechanicalSolution" , "displacement", M_emSolver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
if ( 0 == M_emSolver.comm()->MyPID() )
{
std::cout << "\nLoad from restart: " << restartInput << ", nIterCirculation = " << nIter << ", time = " << t << std::endl;
}
for ( unsigned int i = 0; i < M_emSolver.electroSolverPtr()->globalSolution().size() ; ++i )
{
ElectrophysiologyUtility::importScalarField (M_emSolver.electroSolverPtr()->globalSolution().at(i), "ElectroSolution" , ("Variable" + std::to_string(i)), M_emSolver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
}
if ( 0 == M_emSolver.comm()->MyPID() )
{
std::cout << "\nLoad from restart: " << restartInput << ", nIterCirculation = " << nIter << ", time = " << t << std::endl;
}
ElectrophysiologyUtility::importScalarField (M_emSolver.activationModelPtr() -> fiberActivationPtr(), "ActivationSolution" , "Activation", M_emSolver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
//ElectrophysiologyUtility::importScalarField (solver.activationTimePtr(), "ActivationTimeSolution" , "Activation Time", solver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
if ( 0 == M_emSolver.comm()->MyPID() )
{
std::cout << "\nLoad from restart: " << restartInput << ", nIterCirculation = " << nIter << ", time = " << t << std::endl;
}
circulationSolver.restartFromFile ( restartDir + "solution.dat" , nIter );
}
protected:
......@@ -231,6 +289,17 @@ protected:
std::string 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() ;
};
Real patchForce (const Real& t, const Real& Tmax, const Real& tmax, const Real& tduration) const
{
bool time ( fmod(t-tmax+0.5*tduration, 700.) < tduration && fmod(t-tmax+0.5*tduration, 700.) > 0);
......
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