diff --git a/lecture/mm/heavyoil/sagd/problem.hh b/lecture/mm/heavyoil/sagd/problem.hh index 781517e290d4b45c9ea8e8608702207d0fc39d6e..c87211901dae447682d81c94a013491a9d3b6990 100644 --- a/lecture/mm/heavyoil/sagd/problem.hh +++ b/lecture/mm/heavyoil/sagd/problem.hh @@ -88,6 +88,11 @@ public: totalMassProducedWater_ = 0.0; } + void setEpisodeIdx(Scalar epiIdx ) + { + episodeIdx_ = epiIdx; + } + /*! * \name Boundary conditions */ @@ -160,7 +165,7 @@ public: values[contiWEqIdx] = -0.193;//*0.5; // (55.5 mol*12.5)/3600 mol/s m = 0.193 values[energyEqIdx] = -9132;//*0.5; // J/sec m 9132 } - else if (globalPos[1] > 2.5 - eps_ && globalPos[1] < 3.5 + eps_) // production well + else if (globalPos[1] > (3.1 - eps_) && globalPos[1] < (3.2 + eps_)) // production well { const Scalar elemPressW = elemVolVars[scvf.insideScvIdx()].pressure(wPhaseIdx); //Pressures @@ -184,15 +189,16 @@ public: using std::log; //divided by molarMass() of water to convert from kg/m s to mol/m s - const Scalar qW = (((2*3.1415*0.5*4e-14)/(log(effectiveRadius/wellRadius))) * + Scalar qW = (((2*3.1415*0.5*4e-14)/(log(effectiveRadius/wellRadius))) * densityW * elemMobW * ( elemPressW-pOut))/0.01801528; + //divided by molarMass() of HeavyOil to convert from kg/m s to mol/m s - const Scalar qN = (((2*3.1415*0.5*4e-14)/(log(effectiveRadius/wellRadius))) * + Scalar qN = (((2*3.1415*0.5*4e-14)/(log(effectiveRadius/wellRadius))) * densityN * elemMobN * (elemPressN-pOut))/0.35; Scalar qE; //without cooling: - // qE = qW*0.018*enthW + qN*enthN*0.350; + //qE = qW*0.018*enthW + qN*enthN*0.350; //with cooling: see Diplomarbeit Stefan Roll, Sept. 2015 Scalar wT = elemVolVars[scvf.insideScvIdx()].temperature(); // well temperature @@ -200,6 +206,7 @@ public: { qE = qW*0.018*enthW + qN*enthN*0.350 + (wT-495.)*5000.; // ~3x injected enthalpy std::cout<< "Cooling now! Extracted enthalpy: " << qE << std::endl; + std::cout<< "Well temperature: " << wT << std::endl; } else qE = qW*0.018*enthW + qN*enthN*0.350; @@ -254,7 +261,7 @@ public: const int timeStepIndex = timeLoop.timeStepIndex(); if (timeStepIndex == 0 || - timeStepIndex % 100 == 0 || //after every 1000000 secs + timeStepIndex % 1 == 0 || timeLoop.willBeFinished()) { diff --git a/lecture/mm/heavyoil/sagd/properties.hh b/lecture/mm/heavyoil/sagd/properties.hh index b6f203b27daa579c1771f6430acfe2cb3d4efa5d..f24876dc5e4d5acd931a600d45cca5603f84d003 100644 --- a/lecture/mm/heavyoil/sagd/properties.hh +++ b/lecture/mm/heavyoil/sagd/properties.hh @@ -28,6 +28,7 @@ #include #include +//#include "../myh2oheavyoilwithsimpleh2odensities.hh" #include #include #include diff --git a/lecture/mm/heavyoil/sagd/sagd.cc b/lecture/mm/heavyoil/sagd/sagd.cc index dddf6526519290245ee8bf73e8eff6678cbca590..d1b6da09a7491415e985185ca988d42f080ec65f 100644 --- a/lecture/mm/heavyoil/sagd/sagd.cc +++ b/lecture/mm/heavyoil/sagd/sagd.cc @@ -138,7 +138,9 @@ int main(int argc, char** argv) try vtkWriter.write(0.0); // instantiate time loop + auto episodeIdx = 0.0; auto timeLoop = std::make_shared>(0, dt, tEnd); + problem->setEpisodeIdx(episodeIdx); timeLoop->setMaxTimeStepSize(maxDt); timeLoop->setPeriodicCheckPoint(getParam("TimeLoop.EpisodeLength", std::numeric_limits::max())); @@ -180,11 +182,21 @@ int main(int argc, char** argv) try // report statistics of this time step timeLoop->reportTimeStep(); + if (hasParam("TimeLoop.EpisodeLength")) + if(timeLoop->isCheckPoint()) + { + std::cout << "\n Episode " << episodeIdx << " done \n" << std::endl; + + episodeIdx++; + problem->setEpisodeIdx(episodeIdx); + + // write vtk output + vtkWriter.write(timeLoop->time()); + } + // set new dt as suggested by the newton solver timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize())); - // write vtk output - vtkWriter.write(timeLoop->time()); } while (!timeLoop->finished()); diff --git a/lecture/mm/heavyoil/sagd/sagd.input b/lecture/mm/heavyoil/sagd/sagd.input index 76255c8299593047722c63044cee2973f2cf2369..ad87e1933ea5186f45e9213127af12a0cb78c309 100644 --- a/lecture/mm/heavyoil/sagd/sagd.input +++ b/lecture/mm/heavyoil/sagd/sagd.input @@ -1,7 +1,8 @@ [TimeLoop] DtInitial = 1800 # [s] TEnd = 31467600 # [s] a year -MaxTimeStepSize = 1800 # [s] +MaxTimeStepSize = 86400 # [s] +EpisodeLength = 86400 [Grid] UpperRight = 60 40 @@ -14,8 +15,8 @@ Name = sagd # name passed to the output routines MaxTimeStepDivisions = 100 [Newton] -MaxSteps = 8 -RelTolerance = 1e-02 +MaxRelativeShift = 1e-5 +MaxSteps = 5 [Component] SolidDensity = 2650 diff --git a/lecture/mm/heavyoil/sagdcyclic/problem.hh b/lecture/mm/heavyoil/sagdcyclic/problem.hh index c006c05ddecf510e3c9e90100147ee21c77701f1..b58a99c2429b143808486aad294a64963d5613ea 100644 --- a/lecture/mm/heavyoil/sagdcyclic/problem.hh +++ b/lecture/mm/heavyoil/sagdcyclic/problem.hh @@ -78,6 +78,7 @@ class SagdCyclicProblem : public PorousMediumFlowProblem using FVElementGeometry = typename GetPropType::LocalView; using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; using GlobalPosition = typename SubControlVolumeFace::GlobalPosition; + using TimeLoopPtr = std::shared_ptr>; public: @@ -163,20 +164,24 @@ public: // negative values for injection at injection well if (globalPos[1] > 8.5 - eps_ && globalPos[1] < 9.5 + eps_) { - if (episodeIdx_ % 2 != 0) // if episodeIndex is uneven >> noInjection Phase - { - values[contiNEqIdx] = 0.0; - values[contiWEqIdx] = 0.0; - values[energyEqIdx] = 0.0; - } - else - { - values[contiNEqIdx] = -0.0; - values[contiWEqIdx] = -0.193*2; // (55.5 mol*12.5)/3600 mol/s m = 0.193 - values[energyEqIdx] = -9132*2; // J/sec m 9132 - } + values[contiNEqIdx] = 0.0; + values[contiWEqIdx] = -(0.193*std::sin(7.272205e-5*time()-21600.)+0.193); + values[energyEqIdx] = -(9132.*std::sin(7.272205e-5*time()-21600.)+9132.); + + //if (episodeIdx_ % 2 != 0) // if episodeIndex is uneven >> noInjection Phase + //{ + //values[contiNEqIdx] = 0.0; + //values[contiWEqIdx] = 0.0; + //values[energyEqIdx] = 0.0; + //} + //else + //{ + //values[contiNEqIdx] = -0.0; + //values[contiWEqIdx] = -0.193*2; // (55.5 mol*12.5)/3600 mol/s m = 0.193 + //values[energyEqIdx] = -9132*2; // J/sec m 9132 + //} } - else if (globalPos[1] > 2.5 - eps_ && globalPos[1] < 3.5 + eps_) // production well + else if (globalPos[1] > (3.1 - eps_) && globalPos[1] < (3.2 + eps_)) // production well { const Scalar elemPressW = elemVolVars[scvf.insideScvIdx()].pressure(wPhaseIdx); //Pressures @@ -278,6 +283,16 @@ public: } } + void setTimeLoop(TimeLoopPtr timeLoop) + { + timeLoop_ = timeLoop; + } + + Scalar time() const + { + return timeLoop_->time(); + } + private: // internal method for the initial condition (reused for the // dirichlet conditions!) @@ -299,6 +314,8 @@ private: int episodeIdx_; + TimeLoopPtr timeLoop_; + Scalar totalMassProducedOil_; Scalar totalMassProducedWater_; mutable Scalar massProducedOil_; diff --git a/lecture/mm/heavyoil/sagdcyclic/sagd_cyclic.cc b/lecture/mm/heavyoil/sagdcyclic/sagd_cyclic.cc index da493b428588f2170d7761a4f7a3df1374e6957d..d36dda2daff3794143eded07180b5575fde82f1c 100644 --- a/lecture/mm/heavyoil/sagdcyclic/sagd_cyclic.cc +++ b/lecture/mm/heavyoil/sagdcyclic/sagd_cyclic.cc @@ -112,10 +112,24 @@ int main(int argc, char** argv) try using Problem = GetPropType; auto problem = std::make_shared(fvGridGeometry); + // get some time loop parameters + using Scalar = GetPropType; + const auto tEnd = getParam("TimeLoop.TEnd"); + const auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); + auto dt = getParam("TimeLoop.DtInitial"); + + // instantiate time loop + auto episodeIdx = 0.0; + auto timeLoop = std::make_shared>(0, dt, tEnd); + problem->setEpisodeIdx(episodeIdx); + timeLoop->setMaxTimeStepSize(maxDt); + timeLoop->setPeriodicCheckPoint(getParam("TimeLoop.EpisodeLength", std::numeric_limits::max())); + // the solution vector using SolutionVector = GetPropType; SolutionVector x(fvGridGeometry->numDofs()); problem->applyInitialSolution(x); + problem->setTimeLoop(timeLoop); auto xOld = x; // the grid variables @@ -123,12 +137,6 @@ int main(int argc, char** argv) try auto gridVariables = std::make_shared(problem, fvGridGeometry); gridVariables->init(x); - // get some time loop parameters - using Scalar = GetPropType; - const auto tEnd = getParam("TimeLoop.TEnd"); - const auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); - auto dt = getParam("TimeLoop.DtInitial"); - // intialize the vtk output module using IOFields = GetPropType; VtkOutputModule vtkWriter(*gridVariables, x, problem->name()); @@ -137,13 +145,6 @@ int main(int argc, char** argv) try IOFields::initOutputModule(vtkWriter); //!< Add model specific output fields vtkWriter.write(0.0); - // instantiate time loop - auto episodeIdx = 0.0; - auto timeLoop = std::make_shared>(0, dt, tEnd); - problem->setEpisodeIdx(episodeIdx); - timeLoop->setMaxTimeStepSize(maxDt); - timeLoop->setPeriodicCheckPoint(getParam("TimeLoop.EpisodeLength", std::numeric_limits::max())); - // the assembler with time loop for instationary problem using Assembler = FVAssembler; auto assembler = std::make_shared(problem, fvGridGeometry, gridVariables, timeLoop, xOld); @@ -190,13 +191,16 @@ int main(int argc, char** argv) try episodeIdx++; problem->setEpisodeIdx(episodeIdx); + + // write vtk output + vtkWriter.write(timeLoop->time()); } // set new dt as suggested by the newton solver timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize())); // write vtk output - vtkWriter.write(timeLoop->time()); + //vtkWriter.write(timeLoop->time()); } while (!timeLoop->finished()); diff --git a/lecture/mm/heavyoil/sagdcyclic/sagd_cyclic.input b/lecture/mm/heavyoil/sagdcyclic/sagd_cyclic.input index a854c1e9cdadc9a821595deaf61161e7f1e0e552..29813588f1e2a4b4f865d347fb428b58a7eede20 100644 --- a/lecture/mm/heavyoil/sagdcyclic/sagd_cyclic.input +++ b/lecture/mm/heavyoil/sagdcyclic/sagd_cyclic.input @@ -1,8 +1,8 @@ [TimeLoop] DtInitial = 900 # [s] TEnd = 31467600 # [s] a year -MaxTimeStepSize = 1800 # [s] -EpisodeLength = 43200 # 43200 # [s] (43200s = 12h) +MaxTimeStepSize = 10800 # [s] +EpisodeLength = 21600 # 43200 # [s] (43200s = 12h) [Grid] UpperRight = 60 40 @@ -15,7 +15,11 @@ Name = sagdCyclic # name passed to the output routines MaxTimeStepDivisions = 100 [Newton] -MaxSteps = 8 +#UseLineSearch +#EnableAbsoluteResidualCriterion = true +#MaxAbsoluteResidual = 1e-5 +MaxSteps = 5 +MaxRelativeShift = 1e-5 [Component] SolidDensity = 2650 diff --git a/lecture/mm/heavyoil/sagdcyclichyst/sagd_cyclic_hyst.input b/lecture/mm/heavyoil/sagdcyclichyst/sagd_cyclic_hyst.input index e639a1caf56b322c4be1b4ddb8f0479a958a81c4..508dbfd04f420d3e5d8bc00adf3119ec2e0ec2fa 100644 --- a/lecture/mm/heavyoil/sagdcyclichyst/sagd_cyclic_hyst.input +++ b/lecture/mm/heavyoil/sagdcyclichyst/sagd_cyclic_hyst.input @@ -1,8 +1,8 @@ [TimeLoop] -DtInitial = 100 # [s] +DtInitial = 900 # [s] TEnd = 31467600 # [s] a year -MaxTimeStepSize = 1800 # [s] -EpisodeLength = 43200 # 43200 # [s] (43200s = 12h) +MaxTimeStepSize = 10800 # [s] +EpisodeLength = 21600 # 43200 # [s] (43200s = 12h) [Grid] UpperRight = 60 40 @@ -15,9 +15,8 @@ Name = sagdCyclicHyst # name passed to the output routines MaxTimeStepDivisions = 100 [Newton] -RelTolerance = 1e-02 -MaxSteps = 8 -WriteConvergence = 0 +MaxSteps = 5 +MaxRelativeShift = 1e-5 [Component] SolidDensity = 2650