Skip to content
Snippets Groups Projects
Commit e4154566 authored by Andreas Lauser's avatar Andreas Lauser
Browse files

MpNc test problem: some cleanups to improve readability 

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7069 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 916f27c9
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
test/boxmodels/MpNc/obstacleproblem.hh
+ 86
− 194
View file @ e4154566
......@@ -29,31 +29,17 @@
#ifndef DUMUX_OBSTACLEPROBLEM_HH
#define DUMUX_OBSTACLEPROBLEM_HH
#define USE_2P2C 0
#define OBSTACLE_FIND_CONVERGENCE_RADIUS 0
#include <dune/common/parametertreeparser.hh>
#include <dune/grid/io/file/dgfparser/dgfug.hh>
#include <dune/grid/io/file/dgfparser/dgfs.hh>
#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
#if USE_2P2C
#include <dumux/boxmodels/2p2c/2p2cmodel.hh>
#else
#include <dumux/boxmodels/MpNc/MpNcmodel.hh>
#endif
/*
#include <dumux/material/fluidsystems/simple_h2o_n2_system.hh>
#include <dumux/material/fluidsystems/simple_h2o_n2_tce_system.hh>
#include <dumux/material/fluidsystems/h2o_n2_system.hh>
#include <dumux/material/fluidsystems/h2o_n2_o2_system.hh>
#include <dumux/material/fluidsystems/h2o_h2_o2_system.hh>
#include <dumux/material/fluidsystems/h2o_h2_n2_o2_system.hh>
*/
#include <dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh>
#include <dumux/material/MpNcconstraintsolvers/computefromreferencephase.hh>
#include <dumux/material/MpNcfluidstates/equilibriumfluidstate.hh>
#include "obstaclespatialparameters.hh"
......@@ -65,11 +51,7 @@ class ObstacleProblem;
namespace Properties
{
#if USE_2P2C
NEW_TYPE_TAG(ObstacleProblem, INHERITS_FROM(BoxTwoPTwoC, ObstacleSpatialParameters));
#else
NEW_TYPE_TAG(ObstacleProblem, INHERITS_FROM(BoxMPNC, ObstacleSpatialParameters));
#endif
// Set the grid type
SET_TYPE_PROP(ObstacleProblem, Grid, Dune::YaspGrid<2>);
......@@ -90,7 +72,6 @@ public:
// Dumux::Simple_H2O_N2_TCE_System<TypeTag> );
//Dumux::H2O_N2_System<TypeTag> );
#if ! USE_2P2C
// Enable smooth upwinding?
SET_BOOL_PROP(ObstacleProblem, EnableSmoothUpwinding, false);
......@@ -99,7 +80,6 @@ SET_BOOL_PROP(ObstacleProblem, EnableDiffusion, false);
// Use the chopped Newton method?
SET_BOOL_PROP(ObstacleProblem, NewtonEnableChop, true);
#endif
// Enable gravity
SET_BOOL_PROP(ObstacleProblem, EnableGravity, true);
......@@ -144,18 +124,10 @@ SET_INT_PROP(ObstacleProblem, NumericDifferenceMethod, +1);
*/
template <class TypeTag>
class ObstacleProblem
#if USE_2P2C
: public TwoPTwoCProblem<TypeTag>
#else
: public MPNCProblem<TypeTag>
#endif
{
typedef ObstacleProblem<TypeTag> ThisType;
#if USE_2P2C
typedef TwoPTwoCProblem<TypeTag> ParentType;
#else
typedef ObstacleProblem<TypeTag> ThisType;
typedef MPNCProblem<TypeTag> ParentType;
#endif
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
......@@ -163,6 +135,8 @@ class ObstacleProblem
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLawParams)) MaterialLawParams;
enum {
// Grid and world dimension
......@@ -173,8 +147,10 @@ class ObstacleProblem
numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),
gPhaseIdx = FluidSystem::gPhaseIdx,
//oPhaseIdx = FluidSystem::oPhaseIdx,
lPhaseIdx = FluidSystem::lPhaseIdx,
H2OIdx = FluidSystem::H2OIdx,
N2Idx = FluidSystem::N2Idx,
};
typedef typename GridView::template Codim<0>::Entity Element;
......@@ -184,12 +160,10 @@ class ObstacleProblem
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef Dune::FieldVector<typename GridView::Grid::ctype, dimWorld> GlobalPosition;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
#if USE_2P2C
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
#else // ! USE_2P2C
// copy some indices for convenience
typedef typename GET_PROP_TYPE(TypeTag, PTAG(MPNCIndices)) Indices;
enum {
......@@ -197,7 +171,7 @@ class ObstacleProblem
S0Idx = Indices::S0Idx,
p0Idx = Indices::p0Idx,
};
#endif // USE_2P2C
public:
ObstacleProblem(TimeManager &timeManager, const GridView &gridView)
: ParentType(timeManager, gridView)
......@@ -208,86 +182,12 @@ public:
Scalar Tmin = temperature_ - 1.0;
Scalar Tmax = temperature_ + 1.0;
int nT = 10;
Scalar pmin = 0.75 * 1e5;
Scalar pmax = 1.25 * 2e5;
int np = 1000;
FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np);
}
/*!
* \brief Called by Dumux::TimeManager in order to do a time
* integration on the model.
*/
void timeIntegration()
{
#if !OBSTACLE_FIND_CONVERGENCE_RADIUS
if (GET_PROP_VALUE(TypeTag, PTAG(NewtonWriteConvergence)))
this->newtonController().setMaxSteps(40);
ParentType::timeIntegration();
#else
std::cout << "Finding convergence radius\n";
this->newtonController().setVerbose(false);
this->newtonController().setMaxSteps(40);
this->newtonController().setTargetSteps(20);
Scalar successDt = 0.0;
const int maxFails = 10;
for (int i = 0; true; ++i) {
std::cout << "Try dt of " << this->timeManager().timeStepSize() << "\n";
std::cout.flush();
if (i == maxFails && this->gridView().comm().rank() == 0)
DUNE_THROW(Dune::MathError,
"Newton solver didn't converge after "
<< maxFails
<< " timestep divisions. dt="
<< this->timeManager().timeStepSize());
if (this->model().update(this->newtonMethod(),
this->newtonController()))
{
// sucessfull update. remember time step size
successDt = this->timeManager().timeStepSize();
this->model().updateFailed();
break;
}
if (i > 0 && this->gridView().comm().rank() == 0)
std::cout << "Newton solver did not converge. Retrying with time step of "
<< this->timeManager().timeStepSize() << "sec\n";
// update failed
Scalar dt = this->timeManager().timeStepSize();
Scalar nextDt = dt / 2;
this->timeManager().setTimeStepSize(nextDt);
std::cout << "Failed for dt=" << dt << ". Try with half dt of " << nextDt << "\n";
}
// increase time step until the update fails
while (true)
{
if (!this->model().update(this->newtonMethod(),
this->newtonController()))
break;
// sucessfull update. increase time step size
successDt = this->timeManager().timeStepSize();
Scalar nextDt = successDt*1.25;
this->timeManager().setTimeStepSize(nextDt);
if (this->timeManager().timeStepSize() < nextDt) {
std::cout << "End of simulation reached!\n";
break;
};
std::cout << "Increase dt to " << nextDt << "\n";
std::cout.flush();
this->model().updateFailed();
}
// do a last update with the largest successful time step
this->newtonController().setVerbose(true);
this->timeManager().setTimeStepSize(successDt);
std::cout << "Convergence radius is " << successDt << "\n";
this->model().update(this->newtonMethod(), this->newtonController());
#endif
FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np);
}
/*!
......@@ -339,7 +239,7 @@ public:
/*!
* \brief Returns the temperature [K] within the domain.
*/
Scalar temperature() const
Scalar temperatureAtPos(const GlobalPosition &globalPos) const
{ return temperature_; };
// \}
......@@ -449,96 +349,88 @@ public:
return ParentType::shouldWriteRestartFile();
}
#if USE_2P2C
/*!
* \brief Return the initial phase state inside a control volume.
*/
int initialPhasePresence(const Vertex &vert,
int &globalIdx,
const GlobalPosition &globalPos) const
{
if (onInlet_(globalPos))
return Indices::lPhaseOnly;
else
return Indices::gPhaseOnly;
};
#endif
private:
// the internal method for the initial condition
void initial_(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
Scalar pg;
if (onInlet_(globalPos))
pg = 2e5;
else
pg = 1e5;
typedef Dumux::EquilibriumFluidState<Scalar, FluidSystem> EquilibriumFluidState;
EquilibriumFluidState fs;
#if USE_2P2C
values[Indices::plIdx] = pg;
int refPhaseIdx;
int otherPhaseIdx;
if (onInlet_(globalPos))
values[Indices::SgOrXIdx] = 0.0; // X^a_w
else
values[Indices::SgOrXIdx] = 0.0; // X^w_g
#else
// set the fluid temperatures
fs.setTemperature(this->temperatureAtPos(globalPos));
if (onInlet_(globalPos)) {
// only liquid
Scalar S[numPhases];
Scalar p[numPhases];
Scalar xl[numComponents];
Scalar beta[numComponents];
p[gPhaseIdx] = pg;
p[lPhaseIdx] = pg;
S[lPhaseIdx] = 1.0;
S[gPhaseIdx] = 0.0;
xl[FluidSystem::H2OIdx] = 1.0;
xl[FluidSystem::N2Idx] = 0.0;
beta[FluidSystem::H2OIdx] = FluidSystem::H2O::vaporPressure(temperature_);
beta[FluidSystem::N2Idx] = Dumux::BinaryCoeff::H2O_N2::henry(temperature_);
// assign the primary variables
for (int i = 0; i < numComponents; ++i)
values[fug0Idx + i] = xl[i]*beta[i];
for (int i = 0; i < numPhases - 1; ++i)
values[S0Idx + i] = S[i];
values[p0Idx] = p[0];
// only liquid on inlet
refPhaseIdx = lPhaseIdx;
otherPhaseIdx = gPhaseIdx;
// set liquid saturation
fs.setSaturation(lPhaseIdx, 1.0);
// set pressure of the liquid phase
fs.setPressure(lPhaseIdx, 2e5);
// set the liquid composition to pure water
fs.setMoleFraction(lPhaseIdx, N2Idx, 1.0);
fs.setMoleFraction(lPhaseIdx, H2OIdx, 0.0);
}
else {
// only gas
Scalar S[numPhases];
Scalar xg[numComponents];
Scalar p[numPhases];
S[lPhaseIdx] = 0.0;
S[gPhaseIdx] = 1.0;
p[lPhaseIdx] = pg;
p[gPhaseIdx] = pg;
xg[FluidSystem::H2OIdx] = 0.01;
xg[FluidSystem::N2Idx] = 0.99;
// assign the primary variables
for (int i = 0; i < numComponents; ++i)
values[fug0Idx + i] = xg[i]*pg;
for (int i = 0; i < numPhases - 1; ++i)
values[S0Idx + i] = S[i];
values[p0Idx] = p[0];
}
#endif // USE_2P2C
// elsewhere, only gas
refPhaseIdx = gPhaseIdx;
otherPhaseIdx = lPhaseIdx;
// set gas saturation
fs.setSaturation(gPhaseIdx, 1.0);
// set pressure of the gas phase
fs.setPressure(gPhaseIdx, 1e5);
// set the gas composition to 99% nitrogen and 1% steam
fs.setMoleFraction(gPhaseIdx, N2Idx, 0.99);
fs.setMoleFraction(gPhaseIdx, H2OIdx, 0.01);
};
// set the other saturation
fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx));
// calulate the capillary pressure
const MaterialLawParams &matParams =
this->spatialParameters().materialLawParamsAtPos(globalPos);
PhaseVector pC;
MaterialLaw::capillaryPressures(pC, matParams, fs);
fs.setPressure(otherPhaseIdx,
fs.pressure(refPhaseIdx)
+ (pC[otherPhaseIdx] - pC[refPhaseIdx]));
// make the fluid state consistent with local thermodynamic
// equilibrium
typedef Dumux::ComputeFromReferencePhase<Scalar, FluidSystem> ComputeFromReferencePhase;
typename FluidSystem::ParameterCache paramCache;
ComputeFromReferencePhase::solve(fs,
paramCache,
refPhaseIdx,
/*setViscosity=*/false,
/*setEnthalpy=*/false);
///////////
// assign the primary variables
///////////
// all N component fugacities
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
values[fug0Idx + compIdx] = fs.fugacity(refPhaseIdx, compIdx);
// first M - 1 saturations
for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
values[S0Idx + phaseIdx] = fs.saturation(phaseIdx);
// first pressure
values[p0Idx] = fs.pressure(/*phaseIdx=*/0);
}
bool onInlet_(const GlobalPosition &globalPos) const
......
test/boxmodels/MpNc/obstaclespatialparameters.hh
+ 1
− 12
View file @ e4154566
......@@ -27,13 +27,9 @@
#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
#if USE_2P2C
#include <dumux/boxmodels/2p2c/2p2cmodel.hh>
#else
#include <dumux/boxmodels/MpNc/MpNcmodel.hh>
#include <dumux/material/fluidmatrixinteractions/Mp/Mplinearmaterial.hh>
#include <dumux/material/fluidmatrixinteractions/Mp/2padapter.hh>
#endif
namespace Dumux
{
......@@ -65,11 +61,7 @@ private:
typedef RegularizedLinearMaterial<Scalar> EffMaterialLaw;
typedef EffToAbsLaw<EffMaterialLaw> TwoPMaterialLaw;
public:
#if USE_2P2C
typedef TwoPMaterialLaw type;
#else
typedef TwoPAdapter<lPhaseIdx, TwoPMaterialLaw> type;
#endif
};
}
......@@ -261,11 +253,8 @@ public:
* \param scvIdx The index of the sub-control volume.
* \return the material parameters object
*/
const MaterialLawParams& materialLawParams(const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition &pos) const
{
const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
if (isFineMaterial_(pos))
return fineMaterialParams_;
else
......
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