Commit 35ea4faf authored by Robert Caulk's avatar Robert Caulk

decouple mechanical and partialsat timesteps

parent 1a6b32e2
Pipeline #139909077 passed with stages
in 130 minutes and 59 seconds
......@@ -119,6 +119,7 @@ namespace CGT {
Real meanInitialPorosity;
bool freezeSaturation;
Real permClamp;
Real manualCrackPerm;
#endif
//Handling imposed pressures/fluxes on elements in the form of {point,value} pairs, IPCells contains the cell handles corresponding to point
......
......@@ -74,6 +74,7 @@ namespace CGT {
using _N::surfaceSolidThroatInPore;
using _N::tesselation;
using BaseFlowSolver::manualCrackPerm;
using BaseFlowSolver::bIntrinsicPerm;
using BaseFlowSolver::checkSphereFacetOverlap;
using BaseFlowSolver::clampKValues;
......
......@@ -499,7 +499,7 @@ namespace CGT {
bool ref = Tri.finite_cells_begin()->info().isvisited;
Real meanK = 0, STDEV = 0, meanRadius = 0, meanDistance = 0;
Real infiniteK = 1e10;
//Real infiniteK = 1e10;
for (VCellIterator cellIt = T[currentTes].cellHandles.begin(); cellIt != T[currentTes].cellHandles.end(); cellIt++) {
CellHandle& cell = *cellIt;
......@@ -589,11 +589,8 @@ namespace CGT {
} else {
if (SeM < 0)
cerr << "negative equivalent saturation, linear system will be unstable" << endl;
Real kM = -kFactor
* exp(bIntrinsicPerm
* (avgPoro
- meanInitialPorosity)); //avgPoroOrig)); consider making all perm relative to the mean instaed of the initial...
Real perm = kM * pow(SeM, nUnsatPerm) * area / distance;
Real kM = -kFactor * exp(bIntrinsicPerm * (avgPoro - meanInitialPorosity)); //avgPoroOrig)); consider making all perm relative to the mean instaed of the initial...
Real perm = kM * pow(SeM, nUnsatPerm) * area / distance;
cell->info().kNorm()[j] = (permClamp > 0 and perm > permClamp) ? permClamp : perm;
//cout << "id " << cell->info().id << " perm " << cell->info().kNorm()[j] << " SeM " << SeM << " kM " << kM << " avgPoro " << avgPoro << endl;
}
......@@ -616,7 +613,9 @@ namespace CGT {
}
} else {
cout << "infinite K1!" << endl;
k = infiniteK;
cell->info().kNorm()[j] = manualCrackPerm;
//k = manualCrackPerm; //infiniteK;
} //Will be corrected in the next loop
if (!neighbourCell->info().isGhost)
(neighbourCell->info().kNorm())[Tri.mirror_index(cell, j)] = (cell->info().kNorm())[j];
......
......@@ -285,7 +285,7 @@ class [email protected][email protected] : public PartialEngine
Real averageSlicePressure(Real posY){return solver->averageSlicePressure(posY);}
Real averagePressure(){return solver->averagePressure();}
void emulateAction(){
virtual void emulateAction(){
scene = Omega::instance().getScene().get();
action();}
#ifdef LINSOLV
......
......@@ -24,8 +24,8 @@ PartialSatClayEngine::~PartialSatClayEngine() { }
void PartialSatClayEngine::action()
{
//if (debug) cout << "Entering partialSatEngineAction "<<endl;
//if (elapsedIters > int(partialSatDT/scene->dt) isActivated=true;
//else isActivated=false;
if (partialSatDT!=0) timeStepControl();
if (!isActivated)
return;
timingDeltas->start();
......@@ -116,10 +116,11 @@ void PartialSatClayEngine::action()
timingDeltas->checkpoint("Update_Volumes");
epsVolCumulative += epsVolMax;
retriangulationLastIter++;
if (partialSatDT==0) retriangulationLastIter++;
if (!updateTriangulation)
updateTriangulation = // If not already set true by another function of by the user, check conditions
(defTolerance > 0 && epsVolCumulative > defTolerance) || (meshUpdateInterval > 0 && retriangulationLastIter > meshUpdateInterval);
(defTolerance > 0 && epsVolCumulative > defTolerance) ||
(meshUpdateInterval > 0 && retriangulationLastIter > meshUpdateInterval);
// remesh everytime a bond break occurs (for DFNFlow-JCFPM coupling)
if (breakControlledRemesh)
......@@ -134,7 +135,7 @@ void PartialSatClayEngine::action()
#endif
if (partialSatEngine)
setCellsDSDP(*solver);
solver->gaussSeidel(scene->dt);
solver->gaussSeidel(partialSatDT==0 ? scene->dt : solverDT);
timingDeltas->checkpoint("Factorize + Solve");
if (partialSatEngine) {
//initializeSaturations(*solver);
......@@ -163,8 +164,10 @@ void PartialSatClayEngine::action()
timingDeltas->checkpoint("forces.sync()");
//if (!decoupleForces) computeViscousForces ( *solver );
timingDeltas->checkpoint("viscous forces");
if (!decoupleForces)
applyForces(*solver);
if (!decoupleForces) {
if (partialSatDT!=0) addPermanentForces(*solver);
else applyForces(*solver);
}
timingDeltas->checkpoint("Applying Forces");
if (debug)
cout << "finished computing forces and applying" << endl;
......@@ -260,7 +263,41 @@ void PartialSatClayEngine::action()
/////// Partial Sat Tools /////////
void PartialSatClayEngine::timeStepControl() {
if (((elapsedIters > int(partialSatDT/scene->dt)) and partialSatDT != 0) or first) {
isActivated = true;
retriangulationLastIter+=elapsedIters;
elapsedIters = 0;
if (first) {
collectedDT = scene->dt;
solverDT = scene->dt;
} else {
solverDT = collectedDT;
collectedDT = 0;
}
if (debug) cout << "using flowtime step =" << solverDT << endl;
} else {
if (partialSatDT != 0) {
elapsedIters++;
collectedDT += scene->dt;
isActivated = true;
}
isActivated = emulatingAction ? true : false;
solverDT = scene->dt;
}
}
void PartialSatClayEngine::addPermanentForces(FlowSolver& flow)
{
typedef typename Solver::FiniteVerticesIterator FiniteVerticesIterator;
FiniteVerticesIterator vertices_end = flow.tesselation().Triangulation().finite_vertices_end();
for (FiniteVerticesIterator V_it = flow.tesselation().Triangulation().finite_vertices_begin(); V_it != vertices_end; V_it++) {
scene->forces.setPermForce(V_it->info().id(), makeVector3r(V_it->info().forces));
}
}
void PartialSatClayEngine::blockLowPoroRegions(FlowSolver& flow)
{
......@@ -378,8 +415,9 @@ void PartialSatClayEngine::updatePorosity(FlowSolver& flow)
if (poro > maxPoroClamp)
poro = maxPoroClamp;
if (!freezeSaturation and directlyModifySatFromPoro) { // updatesaturation with respect to volume change
Real dt = partialSatDT==0 ? scene->dt : solverDT;
Real volWater_o
= (cell->info().volume() - cell->info().dv() * scene->dt) * cell->info().porosity * cell->info().saturation;
= (cell->info().volume() - cell->info().dv() * dt) * cell->info().porosity * cell->info().saturation;
cell->info().saturation
= volWater_o / (poro * cell->info().volume()); // update the saturation with respect to new porosity and volume
}
......@@ -1104,6 +1142,7 @@ void PartialSatClayEngine::initSolver(FlowSolver& flow)
flow.meanInitialPorosity = meanInitialPorosity;
flow.freezeSaturation = freezeSaturation;
flow.permClamp = permClamp;
flow.manualCrackPerm = manualCrackPerm;
}
void PartialSatClayEngine::buildTriangulation(Real pZero, Solver& flow)
......@@ -1233,6 +1272,7 @@ void PartialSatClayEngine::initializeVolumes(FlowSolver& flow)
if (flow.fluidBulkModulus > 0 || thermalEngine || iniVoidVolumes) {
cell->info().invVoidVolume() = 1 / (std::abs(cell->info().volume()) - volumeCorrection * flow.volumeSolidPore(cell));
} else if (partialSatEngine) {
if (cell->info().volume() <= 0) cerr << "cell volume zero, bound to be issues" << endl;
cell->info().invVoidVolume() = 1 / std::abs(cell->info().volume());
}
if (!cell->info().isAlpha and !cell->info().isFictious)
......@@ -1246,7 +1286,7 @@ void PartialSatClayEngine::updateVolumes(FlowSolver& flow)
{
if (debug)
cout << "Updating volumes.............." << endl;
Real invDeltaT = 1 / scene->dt;
Real invDeltaT = 1 / (partialSatDT==0 ? scene->dt : solverDT);
epsVolMax = 0;
Real totVol = 0;
Real totDVol = 0;
......
......@@ -158,6 +158,7 @@ public:
void addIncidentParticleIdsToClumpList(CellHandle cell, std::vector<Body::id_t>& clumpIds);
Body::id_t clump(vector<Body::id_t> ids, unsigned int discretization);
void printPorosityToFile(string file);
void addPermanentForces(FlowSolver& flow);
// fracture network
//void crackCellsAbovePoroThreshold(FlowSolver& flow);
......@@ -174,6 +175,7 @@ public:
void removeForceOnVertices(RTriangulation::Facet_circulator& facet, RTriangulation::Finite_edges_iterator& ed_it);
void circulateFacetstoRemoveForces(RTriangulation::Finite_edges_iterator& edge);
void removeForcesOnBrokenBonds();
void timeStepControl();
// void setPositionsBuffer(bool current);
Real leakOffRate;
Real averageAperture;
......@@ -182,12 +184,19 @@ public:
Real crackArea;
Real crackVolume;
Real totalFractureArea;
Real solverDT;
bool emulatingAction;
virtual void initializeVolumes(FlowSolver& flow);
virtual void updateVolumes(FlowSolver& flow);
virtual void buildTriangulation(Real pZero, Solver& flow);
virtual void initSolver(FlowSolver& flow);
virtual void action();
virtual void emulateAction(){
scene = Omega::instance().getScene().get();
emulatingAction=true;
action();
emulatingAction=false;}
void reloadSolver(FlowSolver& flow) { this->initSolver(flow); }
......@@ -285,7 +294,9 @@ public:
((bool,blockIsoCells,true,,"search for cells that might be surrounded by blocked (minerals or cracks) and block them to avoid numerical instabilities."))
((bool,brokenBondsRemoveCapillaryforces,false,,"if true, broken bonds will also remove any capillary forces associated with the area of the crack"))
((bool,directlyModifySatFromPoro,false,,"if true, changes in porosity are used to directly change porosity. Normally, the water retention curve is taking care of this on its own."))
((Real,partialSatDT,0,,"time step used for partial sat engine. The engine will only activate once every partialSatDT/scene->dt steps. Hydromechanical forces estimated and added as persistant forces to particles during non partial sat time steps."))
((Real,partialSatDT,0,,"time step used for partial sat engine. If >0, the engine will only activate once every partialSatDT/scene->dt steps. Hydromechanical forces estimated and added as persistant forces to particles during non partial sat time steps. This value is not exact, see :yref:`PartialSatClayEngine.collectedDT`"))
((int,elapsedIters,0,,"number of mechanical iters since last flow iter."))
((Real,collectedDT,0,,"this is the exact time step that is computed, it enables the stiffness timestep estimate to change dynamically while maintaining an exact match for the flow timestep"))
,/*PartialSatClayEngineT()*/,
solver = shared_ptr<FlowSolver> (new FlowSolver);
,
......
Markdown is supported
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