From 3b19d8e966d2491172436f14b63ae9caa31c06f1 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Thu, 20 Jan 2011 10:03:55 +0000
Subject: [PATCH] replace tabulators by 4 spaces
in some places the indentation might be slightly messed-up by this
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@5042 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
.../msm/1p2cvs2p/external_interface.hh | 69 +-
appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh | 620 +++++++++---------
appl/lecture/msm/1p2cvs2p/lensproblem2p.hh | 43 +-
dumux/boxmodels/1p/1pfluxvariables.hh | 4 +-
.../2p2cni/2p2cniboundaryvariables.hh | 6 +-
dumux/common/pardiso.hh | 18 +-
dumux/common/start.hh | 2 +-
dumux/decoupled/2p2c/2p2cproblem.hh | 4 +-
dumux/decoupled/2p2c/2p2cproperties.hh | 2 +-
dumux/decoupled/2p2c/dec2p2cfluidstate.hh | 104 +--
dumux/decoupled/2p2c/fvpressure2p2c.hh | 310 ++++-----
.../2p2c/fvpressure2p2cmultiphysics.hh | 304 ++++-----
dumux/decoupled/2p2c/fvtransport2p2c.hh | 60 +-
.../2p2c/fvtransport2p2cmultiphysics.hh | 12 +-
dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh | 10 +-
dumux/decoupled/2p2c/variableclass2p2c.hh | 16 +-
dumux/decoupled/common/decoupledproperties.hh | 4 +-
dumux/decoupled/common/impet.hh | 2 +-
test/decoupled/1p/benchmarkresult.hh | 68 +-
test/decoupled/2p/test_transport_problem.hh | 216 +++---
.../2p2c/test_dec2p2c_spatialparams.hh | 2 +-
test/decoupled/2p2c/test_dec2p2cproblem.hh | 4 +-
.../2p2c/test_multiphysics2p2cproblem.hh | 4 +-
23 files changed, 941 insertions(+), 943 deletions(-)
diff --git a/appl/lecture/msm/1p2cvs2p/external_interface.hh b/appl/lecture/msm/1p2cvs2p/external_interface.hh
index 55af6edb09..e612f390f2 100644
--- a/appl/lecture/msm/1p2cvs2p/external_interface.hh
+++ b/appl/lecture/msm/1p2cvs2p/external_interface.hh
@@ -54,7 +54,7 @@ public:
{
using namespace std;
std::cout
- << "-----> ISP: Interface Soil Properties Initialization ...\n";
+ << "-----> ISP: Interface Soil Properties Initialization ...\n";
//ISP input file defenition
ifstream input;
@@ -65,7 +65,7 @@ public:
cout << "\n";
cout << "-----> ISP: Fatal error! - Data read \n";
cout << "-----> ISP: Could not open the input data file: \""
- << isp_filename << "\n";
+ << isp_filename << "\n";
}
//iPCM input file reading
@@ -81,21 +81,21 @@ public:
input >> reader;
ISP_Permeability = atof(reader);
cout << "-----> ISP: Permeability: " << ISP_Permeability
- << "\n";
+ << "\n";
}
if (reader == string("<FinePermeability>"))
{
input >> reader;
ISP_FinePermeability = atof(reader);
cout << "-----> ISP: Fine permeability: "
- << ISP_FinePermeability << "\n";
+ << ISP_FinePermeability << "\n";
}
if (reader == string("<CoarsePermeability>"))
{
input >> reader;
ISP_CoarsePermeability = atof(reader);
cout << "-----> ISP: Coarse permeability: "
- << ISP_CoarsePermeability << "\n";
+ << ISP_CoarsePermeability << "\n";
}
if (reader == string("<Porosity>"))
{
@@ -108,91 +108,91 @@ public:
input >> reader;
ISP_FinePorosity = atof(reader);
cout << "-----> ISP: Fine porosity: " << ISP_FinePorosity
- << "\n";
+ << "\n";
}
if (reader == string("<CoarsePorosity>"))
{
input >> reader;
ISP_CoarsePorosity = atof(reader);
cout << "-----> ISP: Coarse porosity: " << ISP_CoarsePorosity
- << "\n";
+ << "\n";
}
if (reader == string("<LeakageWellPermeability>"))
{
input >> reader;
ISP_LeakageWellPermeability = atof(reader);
cout << "-----> ISP: Leakage Well Permeability: "
- << ISP_LeakageWellPermeability << "\n";
+ << ISP_LeakageWellPermeability << "\n";
}
if (reader == string("<LongitudinalDispersivity>"))
{
input >> reader;
ISP_LongitudinalDispersivity = atof(reader);
cout << "-----> ISP: Longitudinal dispersivity: "
- << ISP_LongitudinalDispersivity << "\n";
+ << ISP_LongitudinalDispersivity << "\n";
}
if (reader == string("<TransverseDispersivity>"))
{
input >> reader;
ISP_TransverseDispersivity = atof(reader);
cout << "-----> ISP: Transverse dispersivity: "
- << ISP_TransverseDispersivity << "\n";
+ << ISP_TransverseDispersivity << "\n";
}
if (reader == string("<FineBrooksCoreyLambda>"))
{
input >> reader;
ISP_FineBrooksCoreyLambda = atof(reader);
cout << "-----> ISP: Brooks-Corey lambda, fine: "
- << ISP_FineBrooksCoreyLambda << "\n";
+ << ISP_FineBrooksCoreyLambda << "\n";
}
if (reader == string("<FineBrooksCoreyEntryPressure>"))
{
input >> reader;
ISP_FineBrooksCoreyEntryPressure = atof(reader);
cout << "-----> ISP: Brooks-Corey entry pressure, fine: "
- << ISP_FineBrooksCoreyEntryPressure << "\n";
+ << ISP_FineBrooksCoreyEntryPressure << "\n";
}
if (reader == string("<CoarseBrooksCoreyLambda>"))
{
input >> reader;
ISP_CoarseBrooksCoreyLambda = atof(reader);
cout << "-----> ISP: Brooks-Corey lambda, coarse: "
- << ISP_CoarseBrooksCoreyLambda << "\n";
+ << ISP_CoarseBrooksCoreyLambda << "\n";
}
if (reader == string("<CoarseBrooksCoreyEntryPressure>"))
{
input >> reader;
ISP_CoarseBrooksCoreyEntryPressure = atof(reader);
cout << "-----> ISP: Brooks-Corey entry pressure, coarse: "
- << ISP_CoarseBrooksCoreyEntryPressure << "\n";
+ << ISP_CoarseBrooksCoreyEntryPressure << "\n";
}
if (reader == string("<FineResidualSaturationWetting>"))
{
input >> reader;
ISP_FineResidualSaturationWetting = atof(reader);
cout << "-----> ISP: Residual saturation wetting phase, fine: "
- << ISP_FineResidualSaturationWetting << "\n";
+ << ISP_FineResidualSaturationWetting << "\n";
}
if (reader == string("<FineResidualSaturationNonWetting>"))
{
input >> reader;
ISP_FineResidualSaturationNonWetting = atof(reader);
cout << "-----> ISP: Residual saturation nonwetting phase, fine: "
- << ISP_FineResidualSaturationNonWetting << "\n";
+ << ISP_FineResidualSaturationNonWetting << "\n";
}
if (reader == string("<CoarseResidualSaturationWetting>"))
{
input >> reader;
ISP_CoarseResidualSaturationWetting = atof(reader);
cout << "-----> ISP: Residual saturation wetting phase, coarse: "
- << ISP_CoarseResidualSaturationWetting << "\n";
+ << ISP_CoarseResidualSaturationWetting << "\n";
}
if (reader == string("<CoarseResidualSaturationNonWetting>"))
{
input >> reader;
ISP_CoarseResidualSaturationNonWetting = atof(reader);
cout << "-----> ISP: Residual saturation nonwetting phase, coarse: "
- << ISP_CoarseResidualSaturationNonWetting << "\n";
+ << ISP_CoarseResidualSaturationNonWetting << "\n";
}
}
input.close();
@@ -216,7 +216,7 @@ public:
{
using namespace std;
std::cout
- << "-----> IFP: Interface Fluid Properties Initialization ...\n";
+ << "-----> IFP: Interface Fluid Properties Initialization ...\n";
//IFP input file defenition
ifstream input;
@@ -227,7 +227,7 @@ public:
cout << "\n";
cout << "-----> IFP: Fatal error! - Data read \n";
cout << "-----> IFP: Could not open the input data file: \""
- << ifp_filename << "\n";
+ << ifp_filename << "\n";
}
//iPCM input file reading
@@ -244,21 +244,21 @@ public:
input >> reader;
IFP_GasDiffCoeff = atof(reader);
cout << "-----> IFP: Gas Diffusion Coefficient: "
- << IFP_GasDiffCoeff << "\n";
+ << IFP_GasDiffCoeff << "\n";
}
if (reader == string("<CO2ResidualSaturation>"))
{
input >> reader;
IFP_CO2ResidSat = atof(reader);
cout << "-----> IFP: Residual Saturation of CO2: "
- << IFP_CO2ResidSat << "\n";
+ << IFP_CO2ResidSat << "\n";
}
if (reader == string("<MolecularDiffusionCoefficient>"))
{
input >> reader;
IFP_MolecularDiffusionCoefficient = atof(reader);
cout << "-----> IFP: Molecular diffusion coefficient: "
- << IFP_MolecularDiffusionCoefficient << "\n";
+ << IFP_MolecularDiffusionCoefficient << "\n";
}
}
input.close();
@@ -288,8 +288,7 @@ public:
//Initialization of IPP Parameters
{
using namespace std;
- std::cout
- << "-----> IPP: Interface Soil Properties Initialization ...\n";
+ std::cout << "-----> IPP: Interface Soil Properties Initialization ...\n";
//IPP input file defenition
ifstream input;
@@ -300,7 +299,7 @@ public:
cout << "\n";
cout << "-----> IPP: Fatal error! - Data read \n";
cout << "-----> IPP: Could not open the input data file: \""
- << ipp_filename << "\n";
+ << ipp_filename << "\n";
}
//iPCM input file reading
@@ -323,63 +322,63 @@ public:
input >> reader;
IPP_InjectionWellRate = atof(reader);
cout << "-----> IPP: Injection Well Rate: "
- << IPP_InjectionWellRate << "\n";
+ << IPP_InjectionWellRate << "\n";
}
if (reader == string("<InjectionWellWindowSize>"))
{
input >> reader;
IPP_InjectionWindowSize = atof(reader);
cout << "-----> IPP: Injection Well Window Size: "
- << IPP_InjectionWindowSize << "\n";
+ << IPP_InjectionWindowSize << "\n";
}
if (reader == string("<UpperPressure>"))
{
input >> reader;
IPP_UpperPressure = atof(reader);
cout << "-----> IPP: Upper pressure: "
- << IPP_UpperPressure << "\n";
+ << IPP_UpperPressure << "\n";
}
if (reader == string("<LowerPressure>"))
{
input >> reader;
IPP_LowerPressure = atof(reader);
cout << "-----> IPP: Lower pressure: "
- << IPP_LowerPressure << "\n";
+ << IPP_LowerPressure << "\n";
}
if (reader == string("<InfiltrationRate>"))
{
input >> reader;
IPP_InfiltrationRate = atof(reader);
cout << "-----> IPP: Infiltration rate: "
- << IPP_InfiltrationRate << "\n";
+ << IPP_InfiltrationRate << "\n";
}
if (reader == string("<MaxTimeStepSize>"))
{
input >> reader;
IPP_MaxTimeStepSize = atof(reader);
cout << "-----> IPP: Maximum time step size: "
- << IPP_MaxTimeStepSize << "\n";
+ << IPP_MaxTimeStepSize << "\n";
}
if (reader == string("<InfiltrationStartTime>"))
{
input >> reader;
IPP_InfiltrationStartTime = atof(reader);
cout << "-----> IPP: Start time of infiltration: "
- << IPP_InfiltrationStartTime << "\n";
+ << IPP_InfiltrationStartTime << "\n";
}
if (reader == string("<InfiltrationEndTime>"))
{
input >> reader;
IPP_InfiltrationEndTime = atof(reader);
cout << "-----> IPP: End time of infiltration: "
- << IPP_InfiltrationEndTime << "\n";
+ << IPP_InfiltrationEndTime << "\n";
}
if (reader == string("<SimulationNumber>"))
{
input >> reader;
IPP_SimulationNumber = atof(reader);
cout << "-----> IPP: Output Name: "
- << IPP_SimulationNumber << "\n";
+ << IPP_SimulationNumber << "\n";
}
}
input.close();
diff --git a/appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh b/appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh
index effb4f05a9..60e666a484 100644
--- a/appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh
+++ b/appl/lecture/msm/1p2cvs2p/lensproblem1p2c.hh
@@ -44,360 +44,360 @@
namespace Dumux
{
- template <class TypeTag>
- class LensProblem;
+template <class TypeTag>
+class LensProblem;
- //////////
- // Specify the properties for the lens problem
- //////////
- namespace Properties
- {
- NEW_TYPE_TAG(LensProblem, INHERITS_FROM(BoxOnePTwoC));
+//////////
+// Specify the properties for the lens problem
+//////////
+namespace Properties
+{
+NEW_TYPE_TAG(LensProblem, INHERITS_FROM(BoxOnePTwoC));
- // Set the grid type
- SET_PROP(LensProblem, Grid)
- {
+// Set the grid type
+SET_PROP(LensProblem, Grid)
+{
#if HAVE_UG
- typedef Dune::UGGrid<2> type;
+ typedef Dune::UGGrid<2> type;
#else
- typedef Dune::YaspGrid<2> type;
- //typedef Dune::SGrid<2, 2> type;
+ typedef Dune::YaspGrid<2> type;
+ //typedef Dune::SGrid<2, 2> type;
#endif
- };
+};
- SET_PROP(LensProblem, LocalFEMSpace)
- {
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
- enum{dim = GridView::dimension};
+SET_PROP(LensProblem, LocalFEMSpace)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ enum{dim = GridView::dimension};
- public:
- typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for cubes
- // typedef Dune::PDELab::P1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for simplices
- };
+public:
+ typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for cubes
+ // typedef Dune::PDELab::P1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for simplices
+};
- // Set the problem property
- SET_PROP(LensProblem, Problem)
- {
- typedef Dumux::LensProblem<TypeTag> type;
- };
+// Set the problem property
+SET_PROP(LensProblem, Problem)
+{
+ typedef Dumux::LensProblem<TypeTag> type;
+};
- // Set fluid configuration
- SET_PROP(LensProblem, FluidSystem)
- {
- typedef Dumux::WaterContaminant<TypeTag> type;
- };
-
- // Set the spatial parameters
- SET_PROP(LensProblem, SpatialParameters)
- {
- typedef Dumux::LensSpatialParameters1p2c<TypeTag> type;
- };
-
- // Enable gravity
- SET_BOOL_PROP(LensProblem, EnableGravity, false);
-
- // Set specific Newton controller
- SET_PROP(LensProblem, NewtonController)
- {
- typedef Dumux::LensOnePTwoCNewtonController<TypeTag> type;
- };
+// Set fluid configuration
+SET_PROP(LensProblem, FluidSystem)
+{
+ typedef Dumux::WaterContaminant<TypeTag> type;
+};
+
+// Set the spatial parameters
+SET_PROP(LensProblem, SpatialParameters)
+{
+ typedef Dumux::LensSpatialParameters1p2c<TypeTag> type;
+};
+
+// Enable gravity
+SET_BOOL_PROP(LensProblem, EnableGravity, false);
+
+// Set specific Newton controller
+SET_PROP(LensProblem, NewtonController)
+{
+ typedef Dumux::LensOnePTwoCNewtonController<TypeTag> type;
+};
+}
+
+/*!
+ * \ingroup TwoPBoxProblems
+ * \brief Soil decontamination problem where a contaminant infiltrates a fully
+ * water saturated medium.
+ *
+ * The whole problem is symmetric.
+ *
+ * The domain is sized 5m times 4m and features a rectangular lens
+ * with low permeablility which spans from (0.8 m , 2 m) to (4 m, 3 m)
+ * and is surrounded by a medium with higher permability.
+ *
+ * On the top and the bottom of the domain Dirichlet boundary conditions
+ * are used, while Neumann conditions apply on the left and right
+ * boundaries.
+ *
+ * The contaminant is injected at the top boundary from 2.25m to 2.75m at a variable rate.
+ *
+ * The Dirichlet boundaries on the top boundary is different to the bottom pressure.
+ *
+ * This problem uses the \ref TwoPBoxModel.
+ *
+ * This problem should typically simulated until \f$t_{\text{end}} =
+ * 30\,000\;s\f$ is reached. A good choice for the initial time step size
+ * is \f$t_{\text{inital}} = 100\;s\f$.
+ *
+ * To run the simulation execute the following line in shell:
+ * <tt>./lens_1p2c 30000 100</tt>
+ */
+template <class TypeTag >
+class LensProblem : public OnePTwoCBoxProblem<TypeTag>
+{
+ typedef LensProblem<TypeTag> ThisType;
+ typedef OnePTwoCBoxProblem<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
+ enum
+ {
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ // indices of the equations
+ contiEqIdx = Indices::contiEqIdx,
+ transEqIdx = Indices::transEqIdx,
+ // indices of the primary variables
+ pressureIdx = Indices::pressureIdx,
+ x1Idx = Indices::x1Idx
+ };
+
+ 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(TimeManager)) TimeManager;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ LensProblem(TimeManager &timeManager,
+ const GridView &gridView,
+ const GlobalPosition &lowerLeft,
+ const GlobalPosition &upperRight,
+ const GlobalPosition &lensLowerLeft,
+ const GlobalPosition &lensUpperRight)
+ : ParentType(timeManager, gridView)
+ {
+ this->spatialParameters().setLensCoords(lensLowerLeft, lensUpperRight);
+
+ bboxMin_[0] = lowerLeft[0];
+ bboxMin_[1] = lowerLeft[1];
+ bboxMax_[0] = upperRight[0];
+ bboxMax_[1] = upperRight[1];
+
+ //load interface-file
+ Dumux::InterfaceProblemProperties interfaceProbProps("interface1p2c.xml");
+
+ upperPressure_ = interfaceProbProps.IPP_UpperPressure;
+ lowerPressure_ = interfaceProbProps.IPP_LowerPressure;
+ infiltrationRate_ = interfaceProbProps.IPP_InfiltrationRate;
+ infiltrationStartTime_= 1.0e-9;//The infiltrations starts always after the first time step!
+ infiltrationEndTime_= interfaceProbProps.IPP_InfiltrationEndTime;
}
/*!
- * \ingroup TwoPBoxProblems
- * \brief Soil decontamination problem where a contaminant infiltrates a fully
- * water saturated medium.
- *
- * The whole problem is symmetric.
- *
- * The domain is sized 5m times 4m and features a rectangular lens
- * with low permeablility which spans from (0.8 m , 2 m) to (4 m, 3 m)
- * and is surrounded by a medium with higher permability.
- *
- * On the top and the bottom of the domain Dirichlet boundary conditions
- * are used, while Neumann conditions apply on the left and right
- * boundaries.
- *
- * The contaminant is injected at the top boundary from 2.25m to 2.75m at a variable rate.
- *
- * The Dirichlet boundaries on the top boundary is different to the bottom pressure.
- *
- * This problem uses the \ref TwoPBoxModel.
- *
- * This problem should typically simulated until \f$t_{\text{end}} =
- * 30\,000\;s\f$ is reached. A good choice for the initial time step size
- * is \f$t_{\text{inital}} = 100\;s\f$.
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
*
- * To run the simulation execute the following line in shell:
- * <tt>./lens_1p2c 30000 100</tt>
+ * This is used as a prefix for files generated by the simulation.
*/
- template <class TypeTag >
- class LensProblem : public OnePTwoCBoxProblem<TypeTag>
+ const char *name() const
{
- typedef LensProblem<TypeTag> ThisType;
- typedef OnePTwoCBoxProblem<TypeTag> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ std::string simName = "lens-1p2c_run";
+ Dumux::InterfaceProblemProperties interfaceProbProps("interface1p2c.xml");
+ Scalar simNum = interfaceProbProps.IPP_SimulationNumber;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
- enum
- {
- // Grid and world dimension
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld,
-
- // indices of the equations
- contiEqIdx = Indices::contiEqIdx,
- transEqIdx = Indices::transEqIdx,
- // indices of the primary variables
- pressureIdx = Indices::pressureIdx,
- x1Idx = Indices::x1Idx
- };
-
- 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(TimeManager)) TimeManager;
-
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::template Codim<dim>::Entity Vertex;
- typedef typename GridView::Intersection Intersection;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
-
- public:
- LensProblem(TimeManager &timeManager,
- const GridView &gridView,
- const GlobalPosition &lowerLeft,
- const GlobalPosition &upperRight,
- const GlobalPosition &lensLowerLeft,
- const GlobalPosition &lensUpperRight)
- : ParentType(timeManager, gridView)
- {
- this->spatialParameters().setLensCoords(lensLowerLeft, lensUpperRight);
+ return (str(boost::format("%s-%02d")
+ %simName%simNum).c_str());
+ }
- bboxMin_[0] = lowerLeft[0];
- bboxMin_[1] = lowerLeft[1];
- bboxMax_[0] = upperRight[0];
- bboxMax_[1] = upperRight[1];
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return 273.15 + 10; // -> 10°C
+ };
- //load interface-file
- Dumux::InterfaceProblemProperties interfaceProbProps("interface1p2c.xml");
+ // \}
- upperPressure_ = interfaceProbProps.IPP_UpperPressure;
- lowerPressure_ = interfaceProbProps.IPP_LowerPressure;
- infiltrationRate_ = interfaceProbProps.IPP_InfiltrationRate;
- infiltrationStartTime_= 1.0e-9;//The infiltrations starts always after the first time step!
- infiltrationEndTime_= interfaceProbProps.IPP_InfiltrationEndTime;
- }
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
- /*!
- * \name Problem parameters
- */
- // \{
-
- /*!
- * \brief The problem name.
- *
- * This is used as a prefix for files generated by the simulation.
- */
- const char *name() const
- {
- std::string simName = "lens-1p2c_run";
- Dumux::InterfaceProblemProperties interfaceProbProps("interface1p2c.xml");
- Scalar simNum = interfaceProbProps.IPP_SimulationNumber;
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ *
+ * \param values The boundary types for the conservation equations
+ * \param vertex The vertex on the boundary for which the
+ * conditions needs to be specified
+ */
+ void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+ {
+ const GlobalPosition globalPos = vertex.geometry().center();
+
+ if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
+ values.setAllDirichlet();
+ else
+ values.setAllNeumann();
+ if (onInlet_(globalPos))
+ values.setNeumann(transEqIdx);
+ if (onLowerBoundary_(globalPos))
+ values.setNeumann(transEqIdx);
+ }
- return (str(boost::format("%s-%02d")
- %simName%simNum).c_str());
- }
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * control volume.
+ *
+ * \param values The dirichlet values for the primary variables
+ * \param vertex The vertex representing the "half volume on the boundary"
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
+ {
+ const GlobalPosition globalPos = vertex.geometry().center();
- /*!
- * \brief Returns the temperature within the domain.
- *
- * This problem assumes a temperature of 10 degrees Celsius.
- */
- Scalar temperature(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
- {
- return 273.15 + 10; // -> 10°C
- };
-
- // \}
-
- /*!
- * \name Boundary conditions
- */
- // \{
-
- /*!
- * \brief Specifies which kind of boundary condition should be
- * used for which equation on a given boundary segment.
- *
- * \param values The boundary types for the conservation equations
- * \param vertex The vertex on the boundary for which the
- * conditions needs to be specified
- */
- void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+ if (onUpperBoundary_(globalPos))
{
- const GlobalPosition globalPos = vertex.geometry().center();
-
- if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
- values.setAllDirichlet();
- else
- values.setAllNeumann();
- if (onInlet_(globalPos))
- values.setNeumann(transEqIdx);
- if (onLowerBoundary_(globalPos))
- values.setNeumann(transEqIdx);
+ values[pressureIdx] = upperPressure_;
+ values[x1Idx] = 0.0;
}
-
- /*!
- * \brief Evaluate the boundary conditions for a dirichlet
- * control volume.
- *
- * \param values The dirichlet values for the primary variables
- * \param vertex The vertex representing the "half volume on the boundary"
- *
- * For this method, the \a values parameter stores primary variables.
- */
- void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
+ else if (onLowerBoundary_(globalPos))
{
- const GlobalPosition globalPos = vertex.geometry().center();
-
- if (onUpperBoundary_(globalPos))
- {
- values[pressureIdx] = upperPressure_;
- values[x1Idx] = 0.0;
- }
- else if (onLowerBoundary_(globalPos))
- {
- values[pressureIdx] = lowerPressure_;
- values[x1Idx] = 0.0;
- }
- else
- values = 0.0;
+ values[pressureIdx] = lowerPressure_;
+ values[x1Idx] = 0.0;
}
+ else
+ values = 0.0;
+ }
- /*!
- * \brief Evaluate the boundary conditions for a neumann
- * boundary segment.
- *
- * For this method, the \a values parameter stores the mass flux
- * in normal direction of each phase. Negative values mean influx.
- */
- void neumann(PrimaryVariables &values,
- const Element &element,
- const FVElementGeometry &fvElemGeom,
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores the mass flux
+ * in normal direction of each phase. Negative values mean influx.
+ */
+ void neumann(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
const Intersection &isIt,
- int scvIdx,
- int boundaryFaceIdx) const
- {
- const GlobalPosition &globalPos
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos
= element.geometry().corner(scvIdx);
- values = 0.0;
+ values = 0.0;
- const Scalar& time = this->timeManager().time();
+ const Scalar& time = this->timeManager().time();
- if (time >= infiltrationStartTime_ && time <= infiltrationEndTime_)
- {
- if (onInlet_(globalPos))
- values[transEqIdx] = -infiltrationRate_;
- }
+ if (time >= infiltrationStartTime_ && time <= infiltrationEndTime_)
+ {
+ if (onInlet_(globalPos))
+ values[transEqIdx] = -infiltrationRate_;
}
- // \}
-
- /*!
- * \name Volume terms
- */
- // \{
-
- /*!
- * \brief Evaluate the source term for all phases within a given
- * sub-control-volume.
- *
- * For this method, the \a values parameter stores the rate mass
- * generated or annihilate per volume unit. Positive values mean
- * that mass is created, negative ones mean that it vanishes.
- */
- void source(PrimaryVariables &values,
+ }
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * For this method, the \a values parameter stores the rate mass
+ * generated or annihilate per volume unit. Positive values mean
+ * that mass is created, negative ones mean that it vanishes.
+ */
+ void source(PrimaryVariables &values,
const Element &element,
const FVElementGeometry &,
int subControlVolumeIdx) const
- {
- values = Scalar(0.0);
- }
+ {
+ values = Scalar(0.0);
+ }
- /*!
- * \brief Evaluate the initial value for a control volume.
- *
- * For this method, the \a values parameter stores primary
- * variables.
- */
- void initial(PrimaryVariables &values,
- const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
- {
- const GlobalPosition &globalPos
- = element.geometry().corner(scvIdx);
-
- // no contaminant, hydrostatic pressure
- const Scalar depth = this->bboxMax()[1] - globalPos[1];
- const Scalar height = this->bboxMax()[1] - this->bboxMin()[1];
-
- values[contiEqIdx] = upperPressure_ - depth/height*(upperPressure_-lowerPressure_);
- //if (globalPos[1]>3.5 && globalPos[1]<3.75 && globalPos[0]>2.25 && globalPos[0]<2.75 )
- //{ values[transEqIdx] = 0.001;}
- //else
- values[transEqIdx] = 0.0;
- }
- // \}
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * For this method, the \a values parameter stores primary
+ * variables.
+ */
+ void initial(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos
+ = element.geometry().corner(scvIdx);
- private:
- bool onLeftBoundary_(const GlobalPosition &globalPos) const
- {
- return globalPos[0] < this->bboxMin()[0] + eps_;
- }
+ // no contaminant, hydrostatic pressure
+ const Scalar depth = this->bboxMax()[1] - globalPos[1];
+ const Scalar height = this->bboxMax()[1] - this->bboxMin()[1];
- bool onRightBoundary_(const GlobalPosition &globalPos) const
- {
- return globalPos[0] > this->bboxMax()[0] - eps_;
- }
+ values[contiEqIdx] = upperPressure_ - depth/height*(upperPressure_-lowerPressure_);
+ //if (globalPos[1]>3.5 && globalPos[1]<3.75 && globalPos[0]>2.25 && globalPos[0]<2.75 )
+ //{ values[transEqIdx] = 0.001;}
+ //else
+ values[transEqIdx] = 0.0;
+ }
+ // \}
- bool onLowerBoundary_(const GlobalPosition &globalPos) const
- {
- return globalPos[1] < this->bboxMin()[1] + eps_;
- }
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[0] < this->bboxMin()[0] + eps_;
+ }
- bool onUpperBoundary_(const GlobalPosition &globalPos) const
- {
- return globalPos[1] > this->bboxMax()[1] - eps_;
- }
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[0] > this->bboxMax()[0] - eps_;
+ }
- bool onInlet_(const GlobalPosition &globalPos) const
- {
- Scalar width = this->bboxMax()[0] - this->bboxMin()[0];
- Scalar lambda = (this->bboxMax()[0] - globalPos[0])/width;
- return onUpperBoundary_(globalPos)
- && (bboxMax_[0] - 0.45*width)/width > lambda
- && lambda > (bboxMax_[0] - 0.55*width)/width;
- }
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[1] < this->bboxMin()[1] + eps_;
+ }
- static const Scalar eps_ = 3e-6;
- GlobalPosition bboxMin_;
- GlobalPosition bboxMax_;
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[1] > this->bboxMax()[1] - eps_;
+ }
- Scalar upperPressure_;
- Scalar lowerPressure_;
- Scalar infiltrationRate_;
- Scalar infiltrationStartTime_;
- Scalar infiltrationEndTime_;
- };
+ bool onInlet_(const GlobalPosition &globalPos) const
+ {
+ Scalar width = this->bboxMax()[0] - this->bboxMin()[0];
+ Scalar lambda = (this->bboxMax()[0] - globalPos[0])/width;
+ return onUpperBoundary_(globalPos)
+ && (bboxMax_[0] - 0.45*width)/width > lambda
+ && lambda > (bboxMax_[0] - 0.55*width)/width;
+ }
+
+ static const Scalar eps_ = 3e-6;
+ GlobalPosition bboxMin_;
+ GlobalPosition bboxMax_;
+
+ Scalar upperPressure_;
+ Scalar lowerPressure_;
+ Scalar infiltrationRate_;
+ Scalar infiltrationStartTime_;
+ Scalar infiltrationEndTime_;
+};
} //end namespace
#endif
diff --git a/appl/lecture/msm/1p2cvs2p/lensproblem2p.hh b/appl/lecture/msm/1p2cvs2p/lensproblem2p.hh
index 06b8677049..25c311a169 100644
--- a/appl/lecture/msm/1p2cvs2p/lensproblem2p.hh
+++ b/appl/lecture/msm/1p2cvs2p/lensproblem2p.hh
@@ -186,23 +186,22 @@ class LensProblem : public TwoPProblem<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
public:
LensProblem(TimeManager &timeManager,
const GridView &gridView,
- const GlobalPosition &lowerLeft,
- const GlobalPosition &upperRight,
+ const GlobalPosition &lowerLeft,
+ const GlobalPosition &upperRight,
const GlobalPosition &lensLowerLeft,
const GlobalPosition &lensUpperRight)
: ParentType(timeManager, gridView)
{
this->spatialParameters().setLensCoords(lensLowerLeft, lensUpperRight);
- bboxMin_[0] = lowerLeft[0];
- bboxMin_[1] = lowerLeft[1];
- bboxMax_[0] = upperRight[0];
- bboxMax_[1] = upperRight[1];
+ bboxMin_[0] = lowerLeft[0];
+ bboxMin_[1] = lowerLeft[1];
+ bboxMax_[0] = upperRight[0];
+ bboxMax_[1] = upperRight[1];
//load interface-file
Dumux::InterfaceProblemProperties interfaceProbProps("interface2p.xml");
@@ -210,8 +209,8 @@ public:
lowerPressure_ = interfaceProbProps.IPP_LowerPressure;
upperPressure_ = interfaceProbProps.IPP_UpperPressure;
infiltrationRate_ = interfaceProbProps.IPP_InfiltrationRate;
- //infiltrationStartTime_= interfaceProbProps.IPP_InfiltrationStartTime;
- infiltrationStartTime_= 1.0e-9;//The infiltrations starts always after the first time step!
+ //infiltrationStartTime_= interfaceProbProps.IPP_InfiltrationStartTime;
+ infiltrationStartTime_= 1.0e-9;//The infiltrations starts always after the first time step!
infiltrationEndTime_= interfaceProbProps.IPP_InfiltrationEndTime;
}
@@ -232,7 +231,7 @@ public:
Scalar simNum = interfaceProbProps.IPP_SimulationNumber;
return (str(boost::format("%s-%02d")
- %simName%simNum).c_str());
+ %simName%simNum).c_str());
}
/*!
@@ -241,8 +240,8 @@ public:
* This problem assumes a temperature of 10 degrees Celsius.
*/
Scalar temperature(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
{
return 273.15 + 10; // -> 10°C
};
@@ -267,7 +266,7 @@ public:
const GlobalPosition globalPos = vertex.geometry().center();
- if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
+ if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
values.setAllDirichlet();
else
values.setAllNeumann();
@@ -300,7 +299,7 @@ public:
values[SnIdx] = 0.0;
}
else
- values = 0.0;
+ values = 0.0;
// Scalar densityW = this->wettingPhase().density();
//
@@ -345,11 +344,11 @@ public:
const Scalar& time = this->timeManager().time();
- if (time >= infiltrationStartTime_ && time <= infiltrationEndTime_)
- {
- if (onInlet_(globalPos))
- values[contiNEqIdx] = -infiltrationRate_; // kg / (m * s)
- }
+ if (time >= infiltrationStartTime_ && time <= infiltrationEndTime_)
+ {
+ if (onInlet_(globalPos))
+ values[contiNEqIdx] = -infiltrationRate_; // kg / (m * s)
+ }
}
// \}
@@ -423,8 +422,8 @@ private:
{
Scalar width = this->bboxMax()[0] - this->bboxMin()[0];
Scalar lambda = (this->bboxMax()[0] - globalPos[0])/width;
- return onUpperBoundary_(globalPos) && (bboxMax_[0]-0.35*width)/width > lambda && lambda > (bboxMax_[0]-0.55*width)/width;
- }
+ return onUpperBoundary_(globalPos) && (bboxMax_[0]-0.35*width)/width > lambda && lambda > (bboxMax_[0]-0.55*width)/width;
+ }
static const Scalar eps_ = 3e-6;
GlobalPosition bboxMin_;
diff --git a/dumux/boxmodels/1p/1pfluxvariables.hh b/dumux/boxmodels/1p/1pfluxvariables.hh
index df22932d0e..c452c20c9e 100644
--- a/dumux/boxmodels/1p/1pfluxvariables.hh
+++ b/dumux/boxmodels/1p/1pfluxvariables.hh
@@ -26,7 +26,7 @@
* the flux of the fluid over a face of a finite volume for the one-phase model.
*
* This means pressure and temperature gradients, phase densities at
- * the integration point, etc.
+ * the integration point, etc.
*/
#ifndef DUMUX_1P_FLUX_VARIABLES_HH
#define DUMUX_1P_FLUX_VARIABLES_HH
@@ -75,7 +75,7 @@ class OnePFluxVariables
typedef Dune::FieldMatrix<Scalar, dim, dim> Tensor;
public:
- /*
+ /*
* \brief The constructor
*
* \param problem The problem
diff --git a/dumux/boxmodels/2p2cni/2p2cniboundaryvariables.hh b/dumux/boxmodels/2p2cni/2p2cniboundaryvariables.hh
index fbbbed3baf..6856d2068c 100644
--- a/dumux/boxmodels/2p2cni/2p2cniboundaryvariables.hh
+++ b/dumux/boxmodels/2p2cni/2p2cniboundaryvariables.hh
@@ -156,7 +156,7 @@ protected:
ScalarGradient tmp(0.0);
// calculate gradients and secondary variables at IPs of the boundary
- for (int idx = 0;
+ for (int idx = 0;
idx < fvElemGeom_.numVertices;
idx++) // loop over adjacent vertices
{
@@ -198,8 +198,8 @@ protected:
for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
{
- densityAtIP_[phaseIdx] += elemDat[idx].density(phaseIdx)*boundaryFace_->shapeValue[idx];
- molarDensityAtIP_[phaseIdx] += elemDat[idx].molarDensity(phaseIdx)*boundaryFace_->shapeValue[idx];
+ densityAtIP_[phaseIdx] += elemDat[idx].density(phaseIdx)*boundaryFace_->shapeValue[idx];
+ molarDensityAtIP_[phaseIdx] += elemDat[idx].molarDensity(phaseIdx)*boundaryFace_->shapeValue[idx];
}
}
diff --git a/dumux/common/pardiso.hh b/dumux/common/pardiso.hh
index abed548e56..49367ccb54 100644
--- a/dumux/common/pardiso.hh
+++ b/dumux/common/pardiso.hh
@@ -319,24 +319,24 @@ public:
initAndFactor(A);
}
- /*! \brief Prepare the preconditioner.
+ /*! \brief Prepare the preconditioner.
- A solver solves a linear operator equation A(x)=b by applying
+ A solver solves a linear operator equation A(x)=b by applying
one or several steps of the preconditioner. The method pre()
is called before the first apply operation.
b and x are right hand side and solution vector of the linear
system respectively. It may. e.g., scale the system, allocate memory or
compute a (I)LU decomposition.
- Note: The ILU decomposition could also be computed in the constructor
+ Note: The ILU decomposition could also be computed in the constructor
or with a separate method of the derived method if several
linear systems with the same matrix are to be solved.
\param x The left hand side of the equation.
\param b The right hand side of the equation.
- */
+ */
virtual void pre (X& x, Y& b) {}
- /*! \brief Apply one step of the preconditioner to the system \f$ A(v)=d \f$.
+ /*! \brief Apply one step of the preconditioner to the system \f$ A(v)=d \f$.
On entry \f$ v=0 \f$ and \f$ d=b-A(x) \f$ (although this might not be
computed in that way. On exit v contains the update, i.e
@@ -345,7 +345,7 @@ public:
the preconditioner.
\param[out] v The update to be computed
\param d The current defect.
- */
+ */
virtual void apply (X& v, const Y& d)
{
#ifdef HAVE_PARDISO
@@ -389,14 +389,14 @@ public:
#endif
}
- /*! \brief Clean up.
+ /*! \brief Clean up.
- This method is called after the last apply call for the
+ This method is called after the last apply call for the
linear system to be solved. Memory may be deallocated safely
here. x is the solution of the linear equation.
\param x The right hand side of the equation.
- */
+ */
virtual void post (X& x)
{
#ifdef HAVE_PARDISO
diff --git a/dumux/common/start.hh b/dumux/common/start.hh
index 76951d236b..2a3e63cd09 100644
--- a/dumux/common/start.hh
+++ b/dumux/common/start.hh
@@ -121,7 +121,7 @@ int startFromDGF(int argc, char **argv)
std::cerr << "WARNING: THE PROGRAM IS STARTED USING MPI, BUT THE GRID IMPLEMENTATION\n"
<< " YOU HAVE CHOSEN IS NOT PARALLEL!\n";
}
- gridPtr.loadBalance();
+ gridPtr.loadBalance();
}
// instantiate and run the concrete problem
diff --git a/dumux/decoupled/2p2c/2p2cproblem.hh b/dumux/decoupled/2p2c/2p2cproblem.hh
index 4918c8aaf0..81884bcc1c 100644
--- a/dumux/decoupled/2p2c/2p2cproblem.hh
+++ b/dumux/decoupled/2p2c/2p2cproblem.hh
@@ -52,8 +52,8 @@ class IMPETProblem2P2C : public IMPETProblem<TypeTag, Implementation>
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
// material properties
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
enum {
diff --git a/dumux/decoupled/2p2c/2p2cproperties.hh b/dumux/decoupled/2p2c/2p2cproperties.hh
index 85941d3520..ed6d4350d9 100644
--- a/dumux/decoupled/2p2c/2p2cproperties.hh
+++ b/dumux/decoupled/2p2c/2p2cproperties.hh
@@ -108,7 +108,7 @@ SET_PROP(DecoupledTwoPTwoC, TwoPIndices)
// set fluid/component information
SET_PROP(DecoupledTwoPTwoC, NumPhases) //!< The number of phases in the 2p model is 2
{
- // the property is created in decoupledproperties.hh
+ // the property is created in decoupledproperties.hh
private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
diff --git a/dumux/decoupled/2p2c/dec2p2cfluidstate.hh b/dumux/decoupled/2p2c/dec2p2cfluidstate.hh
index 3ba3911a74..5f9262be32 100644
--- a/dumux/decoupled/2p2c/dec2p2cfluidstate.hh
+++ b/dumux/decoupled/2p2c/dec2p2cfluidstate.hh
@@ -60,8 +60,8 @@ class DecTwoPTwoCFluidState : public FluidState<typename GET_PROP_TYPE(TypeTag,
};
public:
- enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
- numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),};
+ enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
+ numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),};
public:
/*!
@@ -90,11 +90,11 @@ public:
}
else
phasePressure_[wPhaseIdx] = pw;
- phasePressure_[nPhaseIdx] = pw; // as long as capillary pressure is neglected
+ phasePressure_[nPhaseIdx] = pw; // as long as capillary pressure is neglected
temperature_=temperature;
- //mole equilibrium ratios K for in case wPhase is reference phase
+ //mole equilibrium ratios K for in case wPhase is reference phase
double k1 = FluidSystem::activityCoeff(wPhaseIdx, wCompIdx, temperature_, phasePressure_[nPhaseIdx], *this)
/ phasePressure_[nPhaseIdx]; // = p^wComp_vap / p_n
double k2 = FluidSystem::activityCoeff(wPhaseIdx, nCompIdx, temperature_, phasePressure_[nPhaseIdx], *this)
@@ -102,23 +102,23 @@ public:
// get mole fraction from equilibrium konstants
Scalar xw1 = (1. - k2) / (k1 -k2);
- Scalar xn1 = xw1 * k1;
+ Scalar xn1 = xw1 * k1;
- // transform mole to mass fractions
- massfrac_[wPhaseIdx][wCompIdx] = xw1 * FluidSystem::molarMass(wCompIdx)
+ // transform mole to mass fractions
+ massfrac_[wPhaseIdx][wCompIdx] = xw1 * FluidSystem::molarMass(wCompIdx)
/ ( xw1 * FluidSystem::molarMass(wCompIdx) + (1.-xw1) * FluidSystem::molarMass(nCompIdx) );
- massfrac_[nPhaseIdx][wCompIdx] = xn1 * FluidSystem::molarMass(wCompIdx)
+ massfrac_[nPhaseIdx][wCompIdx] = xn1 * FluidSystem::molarMass(wCompIdx)
/ ( xn1 * FluidSystem::molarMass(wCompIdx) + (1.-xn1) * FluidSystem::molarMass(nCompIdx) );
//mass equilibrium ratios
- equilRatio_[nPhaseIdx][wCompIdx] = massfrac_[nPhaseIdx][wCompIdx] / massfrac_[wPhaseIdx][wCompIdx]; // = Xn1 / Xw1 = K1
- equilRatio_[nPhaseIdx][nCompIdx] = (1.-massfrac_[nPhaseIdx][wCompIdx])/ (1.-massfrac_[wPhaseIdx][wCompIdx]); // =(1.-Xn1) / (1.-Xw1) = K2
- equilRatio_[wPhaseIdx][nCompIdx] = equilRatio_[wPhaseIdx][wCompIdx] = 1.;
+ equilRatio_[nPhaseIdx][wCompIdx] = massfrac_[nPhaseIdx][wCompIdx] / massfrac_[wPhaseIdx][wCompIdx]; // = Xn1 / Xw1 = K1
+ equilRatio_[nPhaseIdx][nCompIdx] = (1.-massfrac_[nPhaseIdx][wCompIdx])/ (1.-massfrac_[wPhaseIdx][wCompIdx]); // =(1.-Xn1) / (1.-Xw1) = K2
+ equilRatio_[wPhaseIdx][nCompIdx] = equilRatio_[wPhaseIdx][wCompIdx] = 1.;
- // phase fraction of nPhase [mass/totalmass]
- nu_[nPhaseIdx] = 0;
+ // phase fraction of nPhase [mass/totalmass]
+ nu_[nPhaseIdx] = 0;
// check if there is enough of component 1 to form a phase
if (Z1 > massfrac_[nPhaseIdx][wCompIdx] && Z1 < massfrac_[wPhaseIdx][wCompIdx])
@@ -129,7 +129,7 @@ public:
massfrac_[nPhaseIdx][wCompIdx] = Z1; // hence, assign complete mass soluted into nPhase
massfrac_[wPhaseIdx][wCompIdx] = 1.;
}
- else // (Z1 >= Xw1) => no nPhase
+ else // (Z1 >= Xw1) => no nPhase
{
nu_[nPhaseIdx] = 0; // no second phase
massfrac_[wPhaseIdx][wCompIdx] = Z1;
@@ -137,16 +137,16 @@ public:
}
// complete array of mass fractions
- massfrac_[wPhaseIdx][nCompIdx] = 1. - massfrac_[wPhaseIdx][wCompIdx];
- massfrac_[nPhaseIdx][nCompIdx] = 1. - massfrac_[nPhaseIdx][wCompIdx];
+ massfrac_[wPhaseIdx][nCompIdx] = 1. - massfrac_[wPhaseIdx][wCompIdx];
+ massfrac_[nPhaseIdx][nCompIdx] = 1. - massfrac_[nPhaseIdx][wCompIdx];
- // complete phase mass fractions
- nu_[wPhaseIdx] = 1. - nu_[nPhaseIdx];
+ // complete phase mass fractions
+ nu_[wPhaseIdx] = 1. - nu_[nPhaseIdx];
-// // let FluidSystem compute partial pressures
-// FluidSystem::computePartialPressures(temperature_, phasePressure_[nPhaseIdx], *this);
+// // let FluidSystem compute partial pressures
+// FluidSystem::computePartialPressures(temperature_, phasePressure_[nPhaseIdx], *this);
- // get densities with correct composition
+ // get densities with correct composition
density_[wPhaseIdx] = FluidSystem::phaseDensity(wPhaseIdx,
temperature,
phasePressure_[wPhaseIdx],
@@ -160,11 +160,11 @@ public:
Sw_ /= (nu_[wPhaseIdx]/density_[wPhaseIdx] + nu_[nPhaseIdx]/density_[nPhaseIdx]);
massConcentration_[wCompIdx] =
- poro * (massfrac_[wPhaseIdx][wCompIdx] * Sw_ * density_[wPhaseIdx]
- + massfrac_[nPhaseIdx][wCompIdx] * (1.-Sw_) * density_[nPhaseIdx]);
+ poro * (massfrac_[wPhaseIdx][wCompIdx] * Sw_ * density_[wPhaseIdx]
+ + massfrac_[nPhaseIdx][wCompIdx] * (1.-Sw_) * density_[nPhaseIdx]);
massConcentration_[nCompIdx] =
- poro * (massfrac_[wPhaseIdx][nCompIdx] * Sw_ * density_[wPhaseIdx]
- + massfrac_[nPhaseIdx][nCompIdx] * (1-Sw_) * density_[nPhaseIdx]);
+ poro * (massfrac_[wPhaseIdx][nCompIdx] * Sw_ * density_[wPhaseIdx]
+ + massfrac_[nPhaseIdx][nCompIdx] * (1-Sw_) * density_[nPhaseIdx]);
}
//! a flash routine for 2p2c if the saturation instead of total concentration is known.
@@ -191,26 +191,26 @@ public:
// assign values
Sw_ = sat;
phasePressure_[wPhaseIdx] = pw;
- phasePressure_[nPhaseIdx] = pw; // as long as capillary pressure is neglected
+ phasePressure_[nPhaseIdx] = pw; // as long as capillary pressure is neglected
temperature_=temperature;
- nu_[nPhaseIdx] = nu_[wPhaseIdx] = NAN; //in contrast to the standard update() method, satflash() does not calculate nu.
+ nu_[nPhaseIdx] = nu_[wPhaseIdx] = NAN; //in contrast to the standard update() method, satflash() does not calculate nu.
- //mole equilibrium ratios K for in case wPhase is reference phase
+ //mole equilibrium ratios K for in case wPhase is reference phase
double k1 = FluidSystem::activityCoeff(wPhaseIdx, wCompIdx, temperature_, phasePressure_[nPhaseIdx], *this)
- / phasePressure_[nPhaseIdx];
+ / phasePressure_[nPhaseIdx];
double k2 = FluidSystem::activityCoeff(wPhaseIdx, nCompIdx, temperature_, phasePressure_[nPhaseIdx], *this)
- / phasePressure_[nPhaseIdx];
+ / phasePressure_[nPhaseIdx];
// get mole fraction from equilibrium konstants
Scalar xw1 = (1. - k2) / (k1 -k2);
- Scalar xn1 = xw1 * k1;
+ Scalar xn1 = xw1 * k1;
- // transform mole to mass fractions
- massfrac_[wPhaseIdx][wCompIdx] = xw1 * FluidSystem::molarMass(wCompIdx)
+ // transform mole to mass fractions
+ massfrac_[wPhaseIdx][wCompIdx] = xw1 * FluidSystem::molarMass(wCompIdx)
/ ( xw1 * FluidSystem::molarMass(wCompIdx) + (1.-xw1) * FluidSystem::molarMass(nCompIdx) );
- massfrac_[nPhaseIdx][wCompIdx] = xn1 * FluidSystem::molarMass(wCompIdx)
+ massfrac_[nPhaseIdx][wCompIdx] = xn1 * FluidSystem::molarMass(wCompIdx)
/ ( xn1 * FluidSystem::molarMass(wCompIdx) + (1.-xn1) * FluidSystem::molarMass(nCompIdx) );
// complete array of mass fractions
@@ -218,11 +218,11 @@ public:
massfrac_[nPhaseIdx][nCompIdx] = 1. - massfrac_[nPhaseIdx][wCompIdx];
//mass equilibrium ratios
- equilRatio_[nPhaseIdx][wCompIdx] = massfrac_[nPhaseIdx][wCompIdx] / massfrac_[wPhaseIdx][wCompIdx]; // = Xn1 / Xw1 = K1
- equilRatio_[nPhaseIdx][nCompIdx] = (1.-massfrac_[nPhaseIdx][wCompIdx])/ (1.-massfrac_[wPhaseIdx][wCompIdx]); // =(1.-Xn1) / (1.-Xw1) = K2
- equilRatio_[wPhaseIdx][nCompIdx] = equilRatio_[wPhaseIdx][wCompIdx] = 1.;
+ equilRatio_[nPhaseIdx][wCompIdx] = massfrac_[nPhaseIdx][wCompIdx] / massfrac_[wPhaseIdx][wCompIdx]; // = Xn1 / Xw1 = K1
+ equilRatio_[nPhaseIdx][nCompIdx] = (1.-massfrac_[nPhaseIdx][wCompIdx])/ (1.-massfrac_[wPhaseIdx][wCompIdx]); // =(1.-Xn1) / (1.-Xw1) = K2
+ equilRatio_[wPhaseIdx][nCompIdx] = equilRatio_[wPhaseIdx][wCompIdx] = 1.;
- // get densities with correct composition
+ // get densities with correct composition
density_[wPhaseIdx] = FluidSystem::phaseDensity(wPhaseIdx,
temperature,
phasePressure_[wPhaseIdx],
@@ -233,11 +233,11 @@ public:
*this);
massConcentration_[wCompIdx] =
- poro * (massfrac_[wPhaseIdx][wCompIdx] * Sw_ * density_[wPhaseIdx]
- + massfrac_[nPhaseIdx][wCompIdx] * (1.-Sw_) * density_[nPhaseIdx]);
+ poro * (massfrac_[wPhaseIdx][wCompIdx] * Sw_ * density_[wPhaseIdx]
+ + massfrac_[nPhaseIdx][wCompIdx] * (1.-Sw_) * density_[nPhaseIdx]);
massConcentration_[nCompIdx] =
- poro * (massfrac_[wPhaseIdx][nCompIdx] * Sw_ * density_[wPhaseIdx]
- + massfrac_[nPhaseIdx][nCompIdx] * (1-Sw_) * density_[nPhaseIdx]);
+ poro * (massfrac_[wPhaseIdx][nCompIdx] * Sw_ * density_[wPhaseIdx]
+ + massfrac_[nPhaseIdx][nCompIdx] * (1-Sw_) * density_[nPhaseIdx]);
}
//@}
/*!
@@ -277,11 +277,11 @@ public:
*/
Scalar moleFrac(int phaseIdx, int compIdx) const
{
- // as the moass fractions are calculated, it is used to determine the mole fractions
- double moleFrac_ = ( massfrac_[phaseIdx][compIdx] / FluidSystem::molarMass(compIdx) ); // = moles of compIdx
- moleFrac_ /= ( massfrac_[phaseIdx][wCompIdx] / FluidSystem::molarMass(wCompIdx)
- + massfrac_[phaseIdx][nCompIdx] / FluidSystem::molarMass(nCompIdx) ); // /= total moles in phase
- return moleFrac_;
+ // as the moass fractions are calculated, it is used to determine the mole fractions
+ double moleFrac_ = ( massfrac_[phaseIdx][compIdx] / FluidSystem::molarMass(compIdx) ); // = moles of compIdx
+ moleFrac_ /= ( massfrac_[phaseIdx][wCompIdx] / FluidSystem::molarMass(wCompIdx)
+ + massfrac_[phaseIdx][nCompIdx] / FluidSystem::molarMass(nCompIdx) ); // /= total moles in phase
+ return moleFrac_;
}
/*!
@@ -316,11 +316,11 @@ public:
*/
Scalar partialPressure(int componentIdx) const
{
- if(componentIdx==nCompIdx)
- return phasePressure_[nPhaseIdx]*moleFrac(nPhaseIdx, nCompIdx);
- if(componentIdx == wCompIdx)
- return phasePressure_[nPhaseIdx]*moleFrac(nPhaseIdx, wCompIdx);
- else
+ if(componentIdx==nCompIdx)
+ return phasePressure_[nPhaseIdx]*moleFrac(nPhaseIdx, nCompIdx);
+ if(componentIdx == wCompIdx)
+ return phasePressure_[nPhaseIdx]*moleFrac(nPhaseIdx, wCompIdx);
+ else
DUNE_THROW(Dune::NotImplemented, "component not found in fluidState!");
}
diff --git a/dumux/decoupled/2p2c/fvpressure2p2c.hh b/dumux/decoupled/2p2c/fvpressure2p2c.hh
index 019e31f42a..18d347f447 100644
--- a/dumux/decoupled/2p2c/fvpressure2p2c.hh
+++ b/dumux/decoupled/2p2c/fvpressure2p2c.hh
@@ -146,14 +146,14 @@ public:
//pressure solution routine: update estimate for secants, assemble, solve.
void pressure(bool solveTwice = true)
{
- //pre-transport to estimate update vector
- Scalar dt_estimate = 0.;
- Dune::dinfo << "secant guess"<< std::endl;
- problem_.transportModel().update(-1, dt_estimate, problem_.variables().updateEstimate(), false);
- //last argument false in update() makes shure that this is estimate and no "real" transport step
+ //pre-transport to estimate update vector
+ Scalar dt_estimate = 0.;
+ Dune::dinfo << "secant guess"<< std::endl;
+ problem_.transportModel().update(-1, dt_estimate, problem_.variables().updateEstimate(), false);
+ //last argument false in update() makes shure that this is estimate and no "real" transport step
problem_.variables().updateEstimate() *= problem_.timeManager().timeStepSize();
- assemble(false); Dune::dinfo << "pressure calculation"<< std::endl;
+ assemble(false); Dune::dinfo << "pressure calculation"<< std::endl;
solve();
return;
@@ -481,17 +481,17 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
}
else
{
- // derivatives of the fluid volume with respect to concentration of components, or pressure
- if (dV_[0][globalIdxI] == 0)
- volumeDerivatives(globalPos, *eIt, dV_[wPhaseIdx][globalIdxI][0], dV_[nPhaseIdx][globalIdxI][0], dv_dp[globalIdxI][0]);
+ // derivatives of the fluid volume with respect to concentration of components, or pressure
+ if (dV_[0][globalIdxI] == 0)
+ volumeDerivatives(globalPos, *eIt, dV_[wPhaseIdx][globalIdxI][0], dV_[nPhaseIdx][globalIdxI][0], dv_dp[globalIdxI][0]);
- source[wPhaseIdx] *= dV_[wPhaseIdx][globalIdxI]; // note: dV_[i][1] = dv_dC1 = dV/dm1
- source[nPhaseIdx] *= dV_[nPhaseIdx][globalIdxI];
+ source[wPhaseIdx] *= dV_[wPhaseIdx][globalIdxI]; // note: dV_[i][1] = dv_dC1 = dV/dm1
+ source[nPhaseIdx] *= dV_[nPhaseIdx][globalIdxI];
}
- f_[globalIdxI] = volume * (source[wPhaseIdx] + source[nPhaseIdx]);
+ f_[globalIdxI] = volume * (source[wPhaseIdx] + source[nPhaseIdx]);
/***********************************/
- // get absolute permeability
+ // get absolute permeability
FieldMatrix permeabilityI(problem_.spatialParameters().intrinsicPermeability(globalPos, *eIt));
// get mobilities and fractional flow factors
@@ -585,11 +585,11 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
if (first) // if we are at the very first iteration we can't calculate phase potentials
{
- // get fractional flow factors in neigbor
+ // get fractional flow factors in neigbor
Scalar fractionalWJ = lambdaWJ / (lambdaWJ+ lambdaNWJ);
Scalar fractionalNWJ = lambdaNWJ / (lambdaWJ+ lambdaNWJ);
- // perform central weighting
+ // perform central weighting
Scalar lambda = (lambdaWI + lambdaWJ) * 0.5 + (lambdaNWI + lambdaNWJ) * 0.5;
entry = fabs(lambda*faceArea*(permeability*unitOuterNormal)/(dist));
@@ -599,9 +599,9 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
}
else
{
- // determine volume derivatives
+ // determine volume derivatives
if (dV_[0][globalIdxJ] == 0)
- volumeDerivatives(globalPosNeighbor, *neighborPointer, dV_[wPhaseIdx][globalIdxJ][0], dV_[nPhaseIdx][globalIdxJ][0], dv_dp[globalIdxJ][0]);
+ volumeDerivatives(globalPosNeighbor, *neighborPointer, dV_[wPhaseIdx][globalIdxJ][0], dV_[nPhaseIdx][globalIdxJ][0], dv_dp[globalIdxJ][0]);
dv_dC1 = (dV_[wPhaseIdx][globalIdxI] + dV_[wPhaseIdx][globalIdxJ]) / 2; // dV/dm1= dV/dC^1
dv_dC2 = (dV_[nPhaseIdx][globalIdxI] + dV_[nPhaseIdx][globalIdxJ]) / 2;
@@ -620,7 +620,7 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
//jochen: central weighting for gravity term
densityW = rhoMeanW; densityNW = rhoMeanNW;
- switch (pressureType) //Markus: hab (unitOuterNormal * distVec)/dist hinzugefuegt
+ switch (pressureType) //Markus: hab (unitOuterNormal * distVec)/dist hinzugefuegt
{
case pw:
{
@@ -640,7 +640,7 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
}
case pglobal:
{
- DUNE_THROW(Dune::NotImplemented, "Global pressure not yet implemented for 2p2c");
+ DUNE_THROW(Dune::NotImplemented, "Global pressure not yet implemented for 2p2c");
// potentialW = (problem_.variables().pressure()[globalIdxI]
// - problem_.variables().pressure()[globalIdxJ] - fMeanNW * (pcI - pcJ)) / dist;
// potentialNW = (problem_.variables().pressure()[globalIdxI]
@@ -652,55 +652,55 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
potentialW += densityW * (unitDistVec * gravity);
potentialNW += densityNW * (unitDistVec * gravity);
- //store potentials for further calculations (velocity, saturation, ...)
- problem_.variables().potentialWetting(globalIdxI, isIndex) = potentialW;
- problem_.variables().potentialNonwetting(globalIdxI, isIndex) = potentialNW;
-
- // initialize convenience shortcuts
- Scalar lambdaW, lambdaN;
- Scalar dV_w(0.), dV_n(0.); // dV_a = \sum_k \rho_a * dv/dC^k * X^k_a
- Scalar gV_w(0.), gV_n(0.); // multipaper eq(3.3) zeile 3 analogon dV_w
-
-
- //do the upwinding of the mobility depending on the phase potentials
- if (potentialW >= 0.)
- {
- dV_w = (dv_dC1 * problem_.variables().wet_X1(globalIdxI) + dv_dC2 * (1. - problem_.variables().wet_X1(globalIdxI)));
- lambdaW = problem_.variables().mobilityWetting(globalIdxI);
- gV_w = (graddv_dC1 * problem_.variables().wet_X1(globalIdxI) + graddv_dC2 * (1. - problem_.variables().wet_X1(globalIdxI)));
- dV_w *= densityWI; gV_w *= densityWI;
- }
- else
- {
- dV_w = (dv_dC1 * problem_.variables().wet_X1(globalIdxJ) + dv_dC2 * (1. - problem_.variables().wet_X1(globalIdxJ)));
- lambdaW = problem_.variables().mobilityWetting(globalIdxJ);
- gV_w = (graddv_dC1 * problem_.variables().wet_X1(globalIdxJ) + graddv_dC2 * (1. - problem_.variables().wet_X1(globalIdxJ)));
- dV_w *= densityWJ; gV_w *= densityWJ;
- }
- if (potentialNW >= 0.)
- {
- dV_n = (dv_dC1 * problem_.variables().nonwet_X1(globalIdxI) + dv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxI)));
- lambdaN = problem_.variables().mobilityNonwetting(globalIdxI);
- gV_n = (graddv_dC1 * problem_.variables().nonwet_X1(globalIdxI) + graddv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxI)));
- dV_n *= densityNWI; gV_n *= densityNWI;
- }
- else
- {
- dV_n = (dv_dC1 * problem_.variables().nonwet_X1(globalIdxJ) + dv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxJ)));
- lambdaN = problem_.variables().mobilityNonwetting(globalIdxJ);
- gV_n = (graddv_dC1 * problem_.variables().nonwet_X1(globalIdxJ) + graddv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxJ)));
- dV_n *= densityNWJ; gV_n *= densityNWJ;
- }
-
- //calculate current matrix entry
+ //store potentials for further calculations (velocity, saturation, ...)
+ problem_.variables().potentialWetting(globalIdxI, isIndex) = potentialW;
+ problem_.variables().potentialNonwetting(globalIdxI, isIndex) = potentialNW;
+
+ // initialize convenience shortcuts
+ Scalar lambdaW, lambdaN;
+ Scalar dV_w(0.), dV_n(0.); // dV_a = \sum_k \rho_a * dv/dC^k * X^k_a
+ Scalar gV_w(0.), gV_n(0.); // multipaper eq(3.3) zeile 3 analogon dV_w
+
+
+ //do the upwinding of the mobility depending on the phase potentials
+ if (potentialW >= 0.)
+ {
+ dV_w = (dv_dC1 * problem_.variables().wet_X1(globalIdxI) + dv_dC2 * (1. - problem_.variables().wet_X1(globalIdxI)));
+ lambdaW = problem_.variables().mobilityWetting(globalIdxI);
+ gV_w = (graddv_dC1 * problem_.variables().wet_X1(globalIdxI) + graddv_dC2 * (1. - problem_.variables().wet_X1(globalIdxI)));
+ dV_w *= densityWI; gV_w *= densityWI;
+ }
+ else
+ {
+ dV_w = (dv_dC1 * problem_.variables().wet_X1(globalIdxJ) + dv_dC2 * (1. - problem_.variables().wet_X1(globalIdxJ)));
+ lambdaW = problem_.variables().mobilityWetting(globalIdxJ);
+ gV_w = (graddv_dC1 * problem_.variables().wet_X1(globalIdxJ) + graddv_dC2 * (1. - problem_.variables().wet_X1(globalIdxJ)));
+ dV_w *= densityWJ; gV_w *= densityWJ;
+ }
+ if (potentialNW >= 0.)
+ {
+ dV_n = (dv_dC1 * problem_.variables().nonwet_X1(globalIdxI) + dv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxI)));
+ lambdaN = problem_.variables().mobilityNonwetting(globalIdxI);
+ gV_n = (graddv_dC1 * problem_.variables().nonwet_X1(globalIdxI) + graddv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxI)));
+ dV_n *= densityNWI; gV_n *= densityNWI;
+ }
+ else
+ {
+ dV_n = (dv_dC1 * problem_.variables().nonwet_X1(globalIdxJ) + dv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxJ)));
+ lambdaN = problem_.variables().mobilityNonwetting(globalIdxJ);
+ gV_n = (graddv_dC1 * problem_.variables().nonwet_X1(globalIdxJ) + graddv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxJ)));
+ dV_n *= densityNWJ; gV_n *= densityNWJ;
+ }
+
+ //calculate current matrix entry
entry = faceArea * (lambdaW * dV_w + lambdaN * dV_n)
- - volume / numberOfFaces * (lambdaW * gV_w + lambdaN * gV_n); // randintegral - gebietsintegral
+ - volume / numberOfFaces * (lambdaW * gV_w + lambdaN * gV_n); // randintegral - gebietsintegral
entry *= fabs((permeability*unitOuterNormal)/(dist));
//calculate right hand side
rightEntry = faceArea * (unitOuterNormal * unitDistVec) * (densityW * lambdaW * dV_w + densityNW * lambdaN * dV_n);
- rightEntry -= volume / numberOfFaces * (densityW * lambdaW * gV_w + densityNW * lambdaN * gV_n);
- rightEntry *= (permeability * gravity); // = multipaper eq(3.3) zeile 2+3 ohne (p-p_k)
+ rightEntry -= volume / numberOfFaces * (densityW * lambdaW * gV_w + densityNW * lambdaN * gV_n);
+ rightEntry *= (permeability * gravity); // = multipaper eq(3.3) zeile 2+3 ohne (p-p_k)
} // end !first
@@ -717,7 +717,7 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
/************* boundary face ************************/
else
{
- // get volume derivatives inside the cell
+ // get volume derivatives inside the cell
dv_dC1 = dV_[wPhaseIdx][globalIdxI];
dv_dC2 = dV_[nPhaseIdx][globalIdxI];
@@ -740,8 +740,8 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
Dune::FieldVector<Scalar, dim> permeability(0);
permeabilityI.mv(unitDistVec, permeability);
- // create a fluid state for the boundary
- FluidState BCfluidState;
+ // create a fluid state for the boundary
+ FluidState BCfluidState;
Scalar temperatureBC = problem_.temperature(globalPosFace, *eIt);
@@ -762,7 +762,7 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
if (first)
{
- Scalar lambda = lambdaWI+lambdaNWI;
+ Scalar lambda = lambdaWI+lambdaNWI;
A_[globalIdxI][globalIdxI] += lambda * faceArea * (permeability * unitOuterNormal) / (dist);
Scalar pressBC = problem_.dirichletPress(globalPosFace, *isIt);
f_[globalIdxI] += lambda * faceArea * pressBC * (permeability * unitOuterNormal) / (dist);
@@ -788,7 +788,7 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
else // nothing declared at boundary
{
Scalar satBound = problem_.variables().saturation()[globalIdxI];
- BCfluidState.satFlash(satBound, pressBC, problem_.spatialParameters().porosity(globalPos, *eIt), temperatureBC);
+ BCfluidState.satFlash(satBound, pressBC, problem_.spatialParameters().porosity(globalPos, *eIt), temperatureBC);
Dune::dwarn << "no boundary saturation/concentration specified on boundary pos " << globalPosFace << std::endl;
}
@@ -891,7 +891,7 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
//do the upwinding of the mobility depending on the phase potentials
Scalar lambdaW, lambdaNW;
- Scalar dV_w, dV_n; // gV_a weglassen, da dV/dc am Rand ortsunabhängig angenommen -> am rand nicht bestimmbar -> nur Randintegral ohne Gebietsintegral
+ Scalar dV_w, dV_n; // gV_a weglassen, da dV/dc am Rand ortsunabhängig angenommen -> am rand nicht bestimmbar -> nur Randintegral ohne Gebietsintegral
if (potentialW >= 0.)
{
@@ -905,7 +905,7 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
{
densityW = densityWBound;
dV_w = (dv_dC1 * BCfluidState.massFrac(wPhaseIdx, wCompIdx)
- + dv_dC2 * BCfluidState.massFrac(wPhaseIdx, nCompIdx));
+ + dv_dC2 * BCfluidState.massFrac(wPhaseIdx, nCompIdx));
dV_w *= densityW;
lambdaW = lambdaWBound;
}
@@ -920,7 +920,7 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
{
densityNW = densityNWBound;
dV_n = (dv_dC1 * BCfluidState.massFrac(nPhaseIdx, wCompIdx)
- + dv_dC2 * BCfluidState.massFrac(nPhaseIdx, nCompIdx));
+ + dv_dC2 * BCfluidState.massFrac(nPhaseIdx, nCompIdx));
dV_n *= densityNW;
lambdaNW = lambdaNWBound;
}
@@ -939,7 +939,7 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
A_[globalIdxI][globalIdxI] += entry;
f_[globalIdxI] += entry * pressBound;
f_[globalIdxI] -= rightEntry * (unitOuterNormal * unitDistVec);
- } //end of if(first) ... else{...
+ } //end of if(first) ... else{...
} // end dirichlet
/**********************************
@@ -947,17 +947,17 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
**********************************/
else
{
- Dune::FieldVector<Scalar,2> J = problem_.neumann(globalPosFace, *isIt);
- if (first)
- {
- J[wPhaseIdx] /= densityWI;
- J[nPhaseIdx] /= densityNWI;
- }
- else
- {
- J[wPhaseIdx] *= dv_dC1;
- J[nPhaseIdx] *= dv_dC2;
- }
+ Dune::FieldVector<Scalar,2> J = problem_.neumann(globalPosFace, *isIt);
+ if (first)
+ {
+ J[wPhaseIdx] /= densityWI;
+ J[nPhaseIdx] /= densityNWI;
+ }
+ else
+ {
+ J[wPhaseIdx] *= dv_dC1;
+ J[nPhaseIdx] *= dv_dC2;
+ }
f_[globalIdxI] -= (J[wPhaseIdx] + J[nPhaseIdx]) * faceArea;
@@ -1049,10 +1049,10 @@ void FVPressure2P2C<TypeTag>::solve()
template<class TypeTag>
void FVPressure2P2C<TypeTag>::initialMaterialLaws(bool compositional)
{
- // initialize the fluid system
+ // initialize the fluid system
FluidState fluidState;
- // iterate through leaf grid an evaluate c0 at cell center
+ // iterate through leaf grid an evaluate c0 at cell center
ElementIterator eItEnd = problem_.gridView().template end<0> ();
for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
{
@@ -1070,34 +1070,34 @@ void FVPressure2P2C<TypeTag>::initialMaterialLaws(bool compositional)
// initial conditions
problem_.variables().capillaryPressure(globalIdx) = 0.;
- Scalar pressW = 0;
- Scalar pressNW = 0;
- Scalar sat_0=0.;
+ Scalar pressW = 0;
+ Scalar pressNW = 0;
+ Scalar sat_0=0.;
BoundaryConditions2p2c::Flags ictype = problem_.initFormulation(globalPos, *eIt); // get type of initial condition
if(!compositional) //means that we do the first approximate guess without compositions
{
- // phase pressures are unknown, so start with an exemplary
- Scalar exemplaryPressure = problem_.referencePressure(globalPos, *eIt);
- pressW = pressNW = problem_.variables().pressure()[globalIdx] = exemplaryPressure;
-
- if (ictype == BoundaryConditions2p2c::saturation) // saturation initial condition
- {
- sat_0 = problem_.initSat(globalPos, *eIt);
- fluidState.satFlash(sat_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
- }
- else if (ictype == BoundaryConditions2p2c::concentration) // concentration initial condition
- {
- Scalar Z1_0 = problem_.initConcentration(globalPos, *eIt);
- fluidState.update(Z1_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
- }
+ // phase pressures are unknown, so start with an exemplary
+ Scalar exemplaryPressure = problem_.referencePressure(globalPos, *eIt);
+ pressW = pressNW = problem_.variables().pressure()[globalIdx] = exemplaryPressure;
+
+ if (ictype == BoundaryConditions2p2c::saturation) // saturation initial condition
+ {
+ sat_0 = problem_.initSat(globalPos, *eIt);
+ fluidState.satFlash(sat_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
+ else if (ictype == BoundaryConditions2p2c::concentration) // concentration initial condition
+ {
+ Scalar Z1_0 = problem_.initConcentration(globalPos, *eIt);
+ fluidState.update(Z1_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
}
- else if(compositional) //means we regard compositional effects since we know an estimate pressure field
+ else if(compositional) //means we regard compositional effects since we know an estimate pressure field
{
- //determine phase pressures from primary pressure variable
- switch (pressureType)
- {
+ //determine phase pressures from primary pressure variable
+ switch (pressureType)
+ {
case pw:
{
pressW = problem_.variables().pressure()[globalIdx];
@@ -1110,18 +1110,18 @@ void FVPressure2P2C<TypeTag>::initialMaterialLaws(bool compositional)
pressNW = problem_.variables().pressure()[globalIdx];
break;
}
- }
-
- if (ictype == Dumux::BoundaryConditions2p2c::saturation) // saturation initial condition
- {
- sat_0 = problem_.initSat(globalPos, *eIt);
- fluidState.satFlash(sat_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
- }
- else if (ictype == Dumux::BoundaryConditions2p2c::concentration) // concentration initial condition
- {
- Scalar Z1_0 = problem_.initConcentration(globalPos, *eIt);
- fluidState.update(Z1_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
- }
+ }
+
+ if (ictype == Dumux::BoundaryConditions2p2c::saturation) // saturation initial condition
+ {
+ sat_0 = problem_.initSat(globalPos, *eIt);
+ fluidState.satFlash(sat_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
+ else if (ictype == Dumux::BoundaryConditions2p2c::concentration) // concentration initial condition
+ {
+ Scalar Z1_0 = problem_.initConcentration(globalPos, *eIt);
+ fluidState.update(Z1_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
}
// initialize densities
@@ -1155,9 +1155,9 @@ void FVPressure2P2C<TypeTag>::initialMaterialLaws(bool compositional)
template<class TypeTag>
void FVPressure2P2C<TypeTag>::updateMaterialLaws()
{
- // this method only completes the variables: = old postprocessupdate()
+ // this method only completes the variables: = old postprocessupdate()
- // instantiate a brandnew fluid state object
+ // instantiate a brandnew fluid state object
FluidState fluidState;
//get timestep for error term
@@ -1232,7 +1232,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLaws()
fluidState.update(Z1, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
/*************************************
- * update variables in variableclass
+ * update variables in variableclass
*************************************/
// initialize saturation
problem_.variables().saturation(globalIdx) = fluidState.saturation(wPhaseIdx);
@@ -1243,9 +1243,9 @@ void FVPressure2P2C<TypeTag>::updateMaterialLaws()
// initialize viscosities
problem_.variables().viscosityWetting(globalIdx)
- = FluidSystem::phaseViscosity(wPhaseIdx, temperature_, pressW, fluidState);
+ = FluidSystem::phaseViscosity(wPhaseIdx, temperature_, pressW, fluidState);
problem_.variables().viscosityNonwetting(globalIdx)
- = FluidSystem::phaseViscosity(nPhaseIdx, temperature_, pressNW, fluidState);
+ = FluidSystem::phaseViscosity(nPhaseIdx, temperature_, pressNW, fluidState);
// initialize mobilities
problem_.variables().mobilityWetting(globalIdx) =
@@ -1261,9 +1261,9 @@ void FVPressure2P2C<TypeTag>::updateMaterialLaws()
// initialize densities
problem_.variables().densityWetting(globalIdx)
- = FluidSystem::phaseDensity(wPhaseIdx, temperature_, pressW, fluidState);
+ = FluidSystem::phaseDensity(wPhaseIdx, temperature_, pressW, fluidState);
problem_.variables().densityNonwetting(globalIdx)
- = FluidSystem::phaseDensity(nPhaseIdx, temperature_, pressNW, fluidState);
+ = FluidSystem::phaseDensity(nPhaseIdx, temperature_, pressNW, fluidState);
problem_.spatialParameters().update(fluidState.saturation(wPhaseIdx), *eIt);
@@ -1274,7 +1274,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLaws()
Scalar massn = problem_.variables().numericalDensity(globalIdx, nPhaseIdx) = sumConc * fluidState.phaseMassFraction(nPhaseIdx);
if ((problem_.variables().densityWetting(globalIdx)*problem_.variables().densityNonwetting(globalIdx)) == 0)
- DUNE_THROW(Dune::MathError, "Decoupled2p2c::postProcessUpdate: try to divide by 0 density");
+ DUNE_THROW(Dune::MathError, "Decoupled2p2c::postProcessUpdate: try to divide by 0 density");
Scalar vol = massw / problem_.variables().densityWetting(globalIdx) + massn / problem_.variables().densityNonwetting(globalIdx);
if (dt != 0)
{
@@ -1311,7 +1311,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLaws()
template<class TypeTag>
void FVPressure2P2C<TypeTag>::volumeDerivatives(GlobalPosition globalPos, ElementPointer ep, Scalar& dv_dC1, Scalar& dv_dC2, Scalar& dv_dp)
{
- // cell index
+ // cell index
int globalIdx = problem_.variables().index(*ep);
// get cell temperature
@@ -1320,24 +1320,24 @@ void FVPressure2P2C<TypeTag>::volumeDerivatives(GlobalPosition globalPos, Elemen
// initialize an Fluid state for the update
FluidState updFluidState;
- /**********************************
- * a) get necessary variables
- **********************************/
- //determine phase pressures from primary pressure variable
+ /**********************************
+ * a) get necessary variables
+ **********************************/
+ //determine phase pressures from primary pressure variable
Scalar pressW=0.;
- switch (pressureType)
- {
- case pw:
- {
- pressW = problem_.variables().pressure()[globalIdx];
- break;
- }
- case pn:
- {
- pressW = problem_.variables().pressure()[globalIdx] - problem_.variables().capillaryPressure(globalIdx);
- break;
- }
- }
+ switch (pressureType)
+ {
+ case pw:
+ {
+ pressW = problem_.variables().pressure()[globalIdx];
+ break;
+ }
+ case pn:
+ {
+ pressW = problem_.variables().pressure()[globalIdx] - problem_.variables().capillaryPressure(globalIdx);
+ break;
+ }
+ }
Scalar v_w = 1. / problem_.variables().densityWetting(globalIdx);
Scalar v_g = 1. / problem_.variables().densityNonwetting(globalIdx);
@@ -1350,18 +1350,18 @@ void FVPressure2P2C<TypeTag>::volumeDerivatives(GlobalPosition globalPos, Elemen
// actual fluid volume
Scalar volalt = (m1+m2) * (nuw1 * v_w + (1-nuw1) * v_g);
- /**********************************
- * b) define increments
- **********************************/
+ /**********************************
+ * b) define increments
+ **********************************/
// increments for numerical derivatives
Scalar inc1 = (fabs(problem_.variables().updateEstimate(globalIdx, wCompIdx)) > 1e-8 / v_w) ? problem_.variables().updateEstimate(globalIdx,wCompIdx) : 1e-8/v_w;
Scalar inc2 =(fabs(problem_.variables().updateEstimate(globalIdx, nCompIdx)) > 1e-8 / v_g) ? problem_.variables().updateEstimate(globalIdx,nCompIdx) : 1e-8 / v_g;
Scalar incp = 1e-2;
- /**********************************
- * c) Secant method for derivatives
- **********************************/
+ /**********************************
+ * c) Secant method for derivatives
+ **********************************/
// numerical derivative of fluid volume with respect to pressure
Scalar p_ = pressW + incp;
@@ -1375,17 +1375,17 @@ void FVPressure2P2C<TypeTag>::volumeDerivatives(GlobalPosition globalPos, Elemen
// numerical derivative of fluid volume with respect to mass of component 1
m1 += inc1;
Z1 = m1 / (m1 + m2);
- updFluidState.update(Z1, pressW, problem_.spatialParameters().porosity(globalPos, *ep), temperature_);
- Scalar satt = updFluidState.saturation(wPhaseIdx);
- Scalar nuw = satt / v_w / (satt/v_w + (1-satt)/v_g);
+ updFluidState.update(Z1, pressW, problem_.spatialParameters().porosity(globalPos, *ep), temperature_);
+ Scalar satt = updFluidState.saturation(wPhaseIdx);
+ Scalar nuw = satt / v_w / (satt/v_w + (1-satt)/v_g);
dv_dC1 = ((m1+m2) * (nuw * v_w + (1-nuw) * v_g) - volalt) /inc1;
m1 -= inc1;
// numerical derivative of fluid volume with respect to mass of component 2
m2 += inc2;
Z1 = m1 / (m1 + m2);
- updFluidState.update(Z1, pressW, problem_.spatialParameters().porosity(globalPos, *ep), temperature_);
- satt = updFluidState.saturation(wPhaseIdx);
+ updFluidState.update(Z1, pressW, problem_.spatialParameters().porosity(globalPos, *ep), temperature_);
+ satt = updFluidState.saturation(wPhaseIdx);
nuw = satt / v_w / (satt/v_w + (1-satt)/v_g);
dv_dC2 = ((m1+m2) * (nuw * v_w + (1-nuw) * v_g) - volalt)/ inc2;
m2 -= inc2;
diff --git a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
index a4c02cb071..03c202fbb2 100644
--- a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
+++ b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
@@ -152,14 +152,14 @@ public:
//pressure solution routine: update estimate for secants, assemble, solve.
void pressure(bool solveTwice = true)
{
- //pre-transport to estimate update vector
- Scalar dt_estimate = 0.;
- Dune::dinfo << "secant guess"<< std::endl;
- problem_.transportModel().update(-1, dt_estimate, problem_.variables().updateEstimate(), false);
- //last argument false in update() makes shure that this is estimate and no "real" transport step
+ //pre-transport to estimate update vector
+ Scalar dt_estimate = 0.;
+ Dune::dinfo << "secant guess"<< std::endl;
+ problem_.transportModel().update(-1, dt_estimate, problem_.variables().updateEstimate(), false);
+ //last argument false in update() makes shure that this is estimate and no "real" transport step
problem_.variables().updateEstimate() *= problem_.timeManager().timeStepSize();
- assemble(false); Dune::dinfo << "pressure calculation"<< std::endl;
+ assemble(false); Dune::dinfo << "pressure calculation"<< std::endl;
solve();
return;
@@ -451,17 +451,17 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
}
else
{
- // derivatives of the fluid volume with respect to concentration of components, or pressure
- if (dV_[0][globalIdxI] == 0)
- volumeDerivatives(globalPos, *eIt, dV_[wPhaseIdx][globalIdxI][0], dV_[nPhaseIdx][globalIdxI][0], dv_dp[globalIdxI][0]);
+ // derivatives of the fluid volume with respect to concentration of components, or pressure
+ if (dV_[0][globalIdxI] == 0)
+ volumeDerivatives(globalPos, *eIt, dV_[wPhaseIdx][globalIdxI][0], dV_[nPhaseIdx][globalIdxI][0], dv_dp[globalIdxI][0]);
- source[wPhaseIdx] *= dV_[wPhaseIdx][globalIdxI]; // note: dV_[i][1] = dv_dC1 = dV/dm1
- source[nPhaseIdx] *= dV_[nPhaseIdx][globalIdxI];
+ source[wPhaseIdx] *= dV_[wPhaseIdx][globalIdxI]; // note: dV_[i][1] = dv_dC1 = dV/dm1
+ source[nPhaseIdx] *= dV_[nPhaseIdx][globalIdxI];
}
- f_[globalIdxI] = volume * (source[wPhaseIdx] + source[nPhaseIdx]);
+ f_[globalIdxI] = volume * (source[wPhaseIdx] + source[nPhaseIdx]);
/********************************************************************/
- // get absolute permeability
+ // get absolute permeability
FieldMatrix permeabilityI(problem_.spatialParameters().intrinsicPermeability(globalPos, *eIt));
// get mobilities and fractional flow factors
@@ -555,11 +555,11 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
if (first) // if we are at the very first iteration we can't calculate phase potentials
{
- // get fractional flow factors in neigbor
+ // get fractional flow factors in neigbor
Scalar fractionalWJ = lambdaWJ / (lambdaWJ+ lambdaNWJ);
Scalar fractionalNWJ = lambdaNWJ / (lambdaWJ+ lambdaNWJ);
- // perform central weighting
+ // perform central weighting
Scalar lambda = (lambdaWI + lambdaWJ) * 0.5 + (lambdaNWI + lambdaNWJ) * 0.5;
entry = fabs(lambda*faceArea*(permeability*unitOuterNormal)/(dist));
@@ -638,7 +638,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
//jochen: central weighting for gravity term
densityW = rhoMeanW; densityNW = rhoMeanNW;
- switch (pressureType) //Markus: hab (unitOuterNormal * distVec)/dist hinzugefuegt
+ switch (pressureType) //Markus: hab (unitOuterNormal * distVec)/dist hinzugefuegt
{
case pw:
{
@@ -658,7 +658,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
}
case pglobal:
{
- DUNE_THROW(Dune::NotImplemented, "Global pressure not yet implemented for 2p2c");
+ DUNE_THROW(Dune::NotImplemented, "Global pressure not yet implemented for 2p2c");
// potentialW = (problem_.variables().pressure()[globalIdxI]
// - problem_.variables().pressure()[globalIdxJ] - fMeanNW * (pcI - pcJ)) / dist;
// potentialNW = (problem_.variables().pressure()[globalIdxI]
@@ -669,62 +669,62 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
potentialW += densityW * (unitDistVec * gravity);
potentialNW += densityNW * (unitDistVec * gravity);
- //store potentials for further calculations (velocity, saturation, ...)
- problem_.variables().potentialWetting(globalIdxI, isIndex) = potentialW;
- problem_.variables().potentialNonwetting(globalIdxI, isIndex) = potentialNW;
+ //store potentials for further calculations (velocity, saturation, ...)
+ problem_.variables().potentialWetting(globalIdxI, isIndex) = potentialW;
+ problem_.variables().potentialNonwetting(globalIdxI, isIndex) = potentialNW;
- // initialize convenience shortcuts
- Scalar lambdaW, lambdaN;
- Scalar dV_w(0.), dV_n(0.); // dV_a = \sum_k \rho_a * dv/dC^k * X^k_a
- Scalar gV_w(0.), gV_n(0.); // multipaper eq(3.3) zeile 3 analogon dV_w
+ // initialize convenience shortcuts
+ Scalar lambdaW, lambdaN;
+ Scalar dV_w(0.), dV_n(0.); // dV_a = \sum_k \rho_a * dv/dC^k * X^k_a
+ Scalar gV_w(0.), gV_n(0.); // multipaper eq(3.3) zeile 3 analogon dV_w
- //do the upwinding of the mobility & density depending on the phase potentials
- if (potentialW >= 0.)
- {
- dV_w = (dv_dC1 * problem_.variables().wet_X1(globalIdxI) + dv_dC2 * (1. - problem_.variables().wet_X1(globalIdxI)));
+ //do the upwinding of the mobility & density depending on the phase potentials
+ if (potentialW >= 0.)
+ {
+ dV_w = (dv_dC1 * problem_.variables().wet_X1(globalIdxI) + dv_dC2 * (1. - problem_.variables().wet_X1(globalIdxI)));
dV_w *= densityWI;
- lambdaW = problem_.variables().mobilityWetting(globalIdxI);
-
- if (problem_.variables().subdomain(globalIdxJ)== 2)
- {
- gV_w = (graddv_dC1 * problem_.variables().wet_X1(globalIdxI) + graddv_dC2 * (1. - problem_.variables().wet_X1(globalIdxI)));
- gV_w *= densityWI;
- }
- }
- else
- {
- dV_w = (dv_dC1 * problem_.variables().wet_X1(globalIdxJ) + dv_dC2 * (1. - problem_.variables().wet_X1(globalIdxJ)));
- dV_w *= densityWJ;
- lambdaW = problem_.variables().mobilityWetting(globalIdxJ);
- if (problem_.variables().subdomain(globalIdxJ)== 2)
- {
- gV_w = (graddv_dC1 * problem_.variables().wet_X1(globalIdxJ) + graddv_dC2 * (1. - problem_.variables().wet_X1(globalIdxJ)));
- gV_w *= densityWJ;
- }
- }
- if (potentialNW >= 0.)
- {
- dV_n = (dv_dC1 * problem_.variables().nonwet_X1(globalIdxI) + dv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxI)));
- dV_n *= densityNWI;
- lambdaN = problem_.variables().mobilityNonwetting(globalIdxI);
- if (problem_.variables().subdomain(globalIdxJ)== 2)
- {
- gV_n = (graddv_dC1 * problem_.variables().nonwet_X1(globalIdxI) + graddv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxI)));
- gV_n *= densityNWI;
- }
- }
- else
- {
- dV_n = (dv_dC1 * problem_.variables().nonwet_X1(globalIdxJ) + dv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxJ)));
- dV_n *= densityNWJ;
- lambdaN = problem_.variables().mobilityNonwetting(globalIdxJ);
- if (problem_.variables().subdomain(globalIdxJ)== 2)
- {
- gV_n = (graddv_dC1 * problem_.variables().nonwet_X1(globalIdxJ) + graddv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxJ)));
- gV_n *= densityNWJ;
- }
- }
+ lambdaW = problem_.variables().mobilityWetting(globalIdxI);
+
+ if (problem_.variables().subdomain(globalIdxJ)== 2)
+ {
+ gV_w = (graddv_dC1 * problem_.variables().wet_X1(globalIdxI) + graddv_dC2 * (1. - problem_.variables().wet_X1(globalIdxI)));
+ gV_w *= densityWI;
+ }
+ }
+ else
+ {
+ dV_w = (dv_dC1 * problem_.variables().wet_X1(globalIdxJ) + dv_dC2 * (1. - problem_.variables().wet_X1(globalIdxJ)));
+ dV_w *= densityWJ;
+ lambdaW = problem_.variables().mobilityWetting(globalIdxJ);
+ if (problem_.variables().subdomain(globalIdxJ)== 2)
+ {
+ gV_w = (graddv_dC1 * problem_.variables().wet_X1(globalIdxJ) + graddv_dC2 * (1. - problem_.variables().wet_X1(globalIdxJ)));
+ gV_w *= densityWJ;
+ }
+ }
+ if (potentialNW >= 0.)
+ {
+ dV_n = (dv_dC1 * problem_.variables().nonwet_X1(globalIdxI) + dv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxI)));
+ dV_n *= densityNWI;
+ lambdaN = problem_.variables().mobilityNonwetting(globalIdxI);
+ if (problem_.variables().subdomain(globalIdxJ)== 2)
+ {
+ gV_n = (graddv_dC1 * problem_.variables().nonwet_X1(globalIdxI) + graddv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxI)));
+ gV_n *= densityNWI;
+ }
+ }
+ else
+ {
+ dV_n = (dv_dC1 * problem_.variables().nonwet_X1(globalIdxJ) + dv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxJ)));
+ dV_n *= densityNWJ;
+ lambdaN = problem_.variables().mobilityNonwetting(globalIdxJ);
+ if (problem_.variables().subdomain(globalIdxJ)== 2)
+ {
+ gV_n = (graddv_dC1 * problem_.variables().nonwet_X1(globalIdxJ) + graddv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxJ)));
+ gV_n *= densityNWJ;
+ }
+ }
//calculate current matrix entry
entry = faceArea * (lambdaW * dV_w + lambdaN * dV_n);
@@ -754,7 +754,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
/************* boundary face ************************/
else
{
- // get volume derivatives inside the cell
+ // get volume derivatives inside the cell
dv_dC1 = dV_[wPhaseIdx][globalIdxI];
dv_dC2 = dV_[nPhaseIdx][globalIdxI];
@@ -777,8 +777,8 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
Dune::FieldVector<Scalar, dim> permeability(0);
permeabilityI.mv(unitDistVec, permeability);
- // create a fluid state for the boundary
- FluidState BCfluidState;
+ // create a fluid state for the boundary
+ FluidState BCfluidState;
Scalar temperatureBC = problem_.temperature(globalPosFace, *eIt);
@@ -799,7 +799,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
if (first)
{
- Scalar lambda = lambdaWI+lambdaNWI;
+ Scalar lambda = lambdaWI+lambdaNWI;
A_[globalIdxI][globalIdxI] += lambda * faceArea * (permeability * unitOuterNormal) / (dist);
pressBC = problem_.dirichletPress(globalPosFace, *isIt);
f_[globalIdxI] += lambda * faceArea * pressBC * (permeability * unitOuterNormal) / (dist);
@@ -825,7 +825,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
else // nothing declared at boundary
{
Scalar satBound = problem_.variables().saturation()[globalIdxI];
- BCfluidState.satFlash(satBound, pressBC, problem_.spatialParameters().porosity(globalPos, *eIt), temperatureBC);
+ BCfluidState.satFlash(satBound, pressBC, problem_.spatialParameters().porosity(globalPos, *eIt), temperatureBC);
Dune::dwarn << "no boundary saturation/concentration specified on boundary pos " << globalPosFace << std::endl;
}
@@ -923,7 +923,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
//do the upwinding of the mobility depending on the phase potentials
Scalar lambdaW, lambdaNW;
- Scalar dV_w, dV_n; // gV_a weglassen, da dV/dc am Rand ortsunabhängig angenommen -> am rand nicht bestimmbar -> nur Randintegral ohne Gebietsintegral
+ Scalar dV_w, dV_n; // gV_a weglassen, da dV/dc am Rand ortsunabhängig angenommen -> am rand nicht bestimmbar -> nur Randintegral ohne Gebietsintegral
if(problem_.variables().subdomain(globalIdxI)==1) // easy 1p subdomain
{
@@ -1003,7 +1003,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
A_[globalIdxI][globalIdxI] += entry;
f_[globalIdxI] += entry * pressBound;
f_[globalIdxI] -= rightEntry * (unitOuterNormal * unitDistVec);
- } //end of if(first) ... else{...
+ } //end of if(first) ... else{...
} // end dirichlet
/**********************************
@@ -1011,17 +1011,17 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
**********************************/
else
{
- Dune::FieldVector<Scalar,2> J = problem_.neumann(globalPosFace, *isIt);
- if (first || problem_.variables().subdomain(globalIdxI)==1)
- {
- J[wPhaseIdx] /= densityWI;
- J[nPhaseIdx] /= densityNWI;
- }
- else
- {
- J[wPhaseIdx] *= dv_dC1;
- J[nPhaseIdx] *= dv_dC2;
- }
+ Dune::FieldVector<Scalar,2> J = problem_.neumann(globalPosFace, *isIt);
+ if (first || problem_.variables().subdomain(globalIdxI)==1)
+ {
+ J[wPhaseIdx] /= densityWI;
+ J[nPhaseIdx] /= densityNWI;
+ }
+ else
+ {
+ J[wPhaseIdx] *= dv_dC1;
+ J[nPhaseIdx] *= dv_dC2;
+ }
f_[globalIdxI] -= (J[wPhaseIdx] + J[nPhaseIdx]) * faceArea;
@@ -1145,10 +1145,10 @@ void FVPressure2P2CMultiPhysics<TypeTag>::solve()
template<class TypeTag>
void FVPressure2P2CMultiPhysics<TypeTag>::initialMaterialLaws(bool compositional)
{
- // initialize the fluid system
+ // initialize the fluid system
FluidState fluidState;
- // iterate through leaf grid an evaluate c0 at cell center
+ // iterate through leaf grid an evaluate c0 at cell center
ElementIterator eItEnd = problem_.gridView().template end<0> ();
for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
{
@@ -1166,34 +1166,34 @@ void FVPressure2P2CMultiPhysics<TypeTag>::initialMaterialLaws(bool compositional
// initial conditions
problem_.variables().capillaryPressure(globalIdx) = 0.;
- Scalar pressW = 0;
- Scalar pressNW = 0;
- Scalar sat_0=0.;
+ Scalar pressW = 0;
+ Scalar pressNW = 0;
+ Scalar sat_0=0.;
BoundaryConditions2p2c::Flags ictype = problem_.initFormulation(globalPos, *eIt); // get type of initial condition
if(!compositional) //means that we do the first approximate guess without compositions
{
- // phase pressures are unknown, so start with an exemplary
- Scalar exemplaryPressure = problem_.referencePressure(globalPos, *eIt);
- pressW = pressNW = problem_.variables().pressure()[globalIdx] = exemplaryPressure;
+ // phase pressures are unknown, so start with an exemplary
+ Scalar exemplaryPressure = problem_.referencePressure(globalPos, *eIt);
+ pressW = pressNW = problem_.variables().pressure()[globalIdx] = exemplaryPressure;
- if (ictype == BoundaryConditions2p2c::saturation) // saturation initial condition
- {
- sat_0 = problem_.initSat(globalPos, *eIt);
- fluidState.satFlash(sat_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
- }
- else if (ictype == BoundaryConditions2p2c::concentration) // concentration initial condition
- {
- Scalar Z1_0 = problem_.initConcentration(globalPos, *eIt);
- fluidState.update(Z1_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
- }
+ if (ictype == BoundaryConditions2p2c::saturation) // saturation initial condition
+ {
+ sat_0 = problem_.initSat(globalPos, *eIt);
+ fluidState.satFlash(sat_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
+ else if (ictype == BoundaryConditions2p2c::concentration) // concentration initial condition
+ {
+ Scalar Z1_0 = problem_.initConcentration(globalPos, *eIt);
+ fluidState.update(Z1_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
}
- else if(compositional) //means we regard compositional effects since we know an estimate pressure field
+ else if(compositional) //means we regard compositional effects since we know an estimate pressure field
{
- //determine phase pressures from primary pressure variable
- switch (pressureType)
- {
+ //determine phase pressures from primary pressure variable
+ switch (pressureType)
+ {
case pw:
{
pressW = problem_.variables().pressure()[globalIdx];
@@ -1206,18 +1206,18 @@ void FVPressure2P2CMultiPhysics<TypeTag>::initialMaterialLaws(bool compositional
pressNW = problem_.variables().pressure()[globalIdx];
break;
}
- }
-
- if (ictype == Dumux::BoundaryConditions2p2c::saturation) // saturation initial condition
- {
- sat_0 = problem_.initSat(globalPos, *eIt);
- fluidState.satFlash(sat_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
- }
- else if (ictype == Dumux::BoundaryConditions2p2c::concentration) // concentration initial condition
- {
- Scalar Z1_0 = problem_.initConcentration(globalPos, *eIt);
- fluidState.update(Z1_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
- }
+ }
+
+ if (ictype == Dumux::BoundaryConditions2p2c::saturation) // saturation initial condition
+ {
+ sat_0 = problem_.initSat(globalPos, *eIt);
+ fluidState.satFlash(sat_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
+ else if (ictype == Dumux::BoundaryConditions2p2c::concentration) // concentration initial condition
+ {
+ Scalar Z1_0 = problem_.initConcentration(globalPos, *eIt);
+ fluidState.update(Z1_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
}
// initialize densities
@@ -1256,9 +1256,9 @@ void FVPressure2P2CMultiPhysics<TypeTag>::initialMaterialLaws(bool compositional
template<class TypeTag>
void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
{
- // this method only completes the variables: = old postprocessupdate()
+ // this method only completes the variables: = old postprocessupdate()
- // instantiate standard 2p2c and pseudo1p fluid state objects
+ // instantiate standard 2p2c and pseudo1p fluid state objects
FluidState fluidState;
PseudoOnePTwoCFluidState<TypeTag> pseudoFluidState;
@@ -1484,7 +1484,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
template<class TypeTag>
void FVPressure2P2CMultiPhysics<TypeTag>::volumeDerivatives(GlobalPosition globalPos, ElementPointer ep, Scalar& dv_dC1, Scalar& dv_dC2, Scalar& dv_dp)
{
- // cell index
+ // cell index
int globalIdx = problem_.variables().index(*ep);
// get cell temperature
@@ -1493,24 +1493,24 @@ void FVPressure2P2CMultiPhysics<TypeTag>::volumeDerivatives(GlobalPosition globa
// initialize an Fluid state for the update
FluidState updFluidState;
- /**********************************
- * a) get necessary variables
- **********************************/
- //determine phase pressures from primary pressure variable
+ /**********************************
+ * a) get necessary variables
+ **********************************/
+ //determine phase pressures from primary pressure variable
Scalar pressW=0.;
- switch (pressureType)
- {
- case pw:
- {
- pressW = problem_.variables().pressure()[globalIdx];
- break;
- }
- case pn:
- {
- pressW = problem_.variables().pressure()[globalIdx] - problem_.variables().capillaryPressure(globalIdx);
- break;
- }
- }
+ switch (pressureType)
+ {
+ case pw:
+ {
+ pressW = problem_.variables().pressure()[globalIdx];
+ break;
+ }
+ case pn:
+ {
+ pressW = problem_.variables().pressure()[globalIdx] - problem_.variables().capillaryPressure(globalIdx);
+ break;
+ }
+ }
Scalar v_w = 1. / problem_.variables().densityWetting(globalIdx);
Scalar v_g = 1. / problem_.variables().densityNonwetting(globalIdx);
@@ -1523,18 +1523,18 @@ void FVPressure2P2CMultiPhysics<TypeTag>::volumeDerivatives(GlobalPosition globa
// actual fluid volume
Scalar volalt = (m1+m2) * (nuw1 * v_w + (1-nuw1) * v_g);
- /**********************************
- * b) define increments
- **********************************/
+ /**********************************
+ * b) define increments
+ **********************************/
// increments for numerical derivatives
Scalar inc1 = (fabs(problem_.variables().updateEstimate(globalIdx, wCompIdx)) > 1e-8 / v_w) ? problem_.variables().updateEstimate(globalIdx,wCompIdx) : 1e-8/v_w;
Scalar inc2 =(fabs(problem_.variables().updateEstimate(globalIdx, nCompIdx)) > 1e-8 / v_g) ? problem_.variables().updateEstimate(globalIdx,nCompIdx) : 1e-8 / v_g;
Scalar incp = 1e-2;
- /**********************************
- * c) Secant method for derivatives
- **********************************/
+ /**********************************
+ * c) Secant method for derivatives
+ **********************************/
// numerical derivative of fluid volume with respect to pressure
Scalar p_ = pressW + incp;
@@ -1548,17 +1548,17 @@ void FVPressure2P2CMultiPhysics<TypeTag>::volumeDerivatives(GlobalPosition globa
// numerical derivative of fluid volume with respect to mass of component 1
m1 += inc1;
Z1 = m1 / (m1 + m2);
- updFluidState.update(Z1, pressW, problem_.spatialParameters().porosity(globalPos, *ep), temperature_);
- Scalar satt = updFluidState.saturation(wPhaseIdx);
- Scalar nuw = satt / v_w / (satt/v_w + (1-satt)/v_g);
+ updFluidState.update(Z1, pressW, problem_.spatialParameters().porosity(globalPos, *ep), temperature_);
+ Scalar satt = updFluidState.saturation(wPhaseIdx);
+ Scalar nuw = satt / v_w / (satt/v_w + (1-satt)/v_g);
dv_dC1 = ((m1+m2) * (nuw * v_w + (1-nuw) * v_g) - volalt) /inc1;
m1 -= inc1;
// numerical derivative of fluid volume with respect to mass of component 2
m2 += inc2;
Z1 = m1 / (m1 + m2);
- updFluidState.update(Z1, pressW, problem_.spatialParameters().porosity(globalPos, *ep), temperature_);
- satt = updFluidState.saturation(wPhaseIdx);
+ updFluidState.update(Z1, pressW, problem_.spatialParameters().porosity(globalPos, *ep), temperature_);
+ satt = updFluidState.saturation(wPhaseIdx);
nuw = satt / v_w / (satt/v_w + (1-satt)/v_g);
dv_dC2 = ((m1+m2) * (nuw * v_w + (1-nuw) * v_g) - volalt)/ inc2;
m2 -= inc2;
diff --git a/dumux/decoupled/2p2c/fvtransport2p2c.hh b/dumux/decoupled/2p2c/fvtransport2p2c.hh
index fe36818af2..17ef95f1ac 100644
--- a/dumux/decoupled/2p2c/fvtransport2p2c.hh
+++ b/dumux/decoupled/2p2c/fvtransport2p2c.hh
@@ -330,8 +330,8 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
// compute mean permeability
Dune::FieldMatrix<Scalar,dim,dim> meanK_(0.);
Dumux::harmonicMeanMatrix(meanK_,
- K_I,
- problem_.spatialParameters().intrinsicPermeability(globalPosNeighbor, *neighborPointer));
+ K_I,
+ problem_.spatialParameters().intrinsicPermeability(globalPosNeighbor, *neighborPointer));
Dune::FieldVector<Scalar,dim> K(0);
meanK_.umv(unitDistVec,K);
@@ -383,7 +383,7 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
}
/******************************************
- * Boundary Face
+ * Boundary Face
******************************************/
if (isIt->boundary())
{
@@ -433,29 +433,29 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
{
Scalar satBound = problem_.dirichletTransport(globalPosFace, *isIt);
BCfluidState.satFlash(satBound, pressBound,
- problem_.spatialParameters().porosity(globalPos, *eIt),
- problem_.temperature(globalPosFace, *eIt));
+ problem_.spatialParameters().porosity(globalPos, *eIt),
+ problem_.temperature(globalPosFace, *eIt));
}
if (bctype == BoundaryConditions2p2c::concentration)
{
// saturation and hence pc and hence corresponding pressure unknown
- Scalar Z1Bound = problem_.dirichletTransport(globalPosFace, *isIt);
- BCfluidState.update(Z1Bound, pressBound,
- problem_.spatialParameters().porosity(globalPos, *eIt),
- problem_.temperature(globalPosFace, *eIt));
+ Scalar Z1Bound = problem_.dirichletTransport(globalPosFace, *isIt);
+ BCfluidState.update(Z1Bound, pressBound,
+ problem_.spatialParameters().porosity(globalPos, *eIt),
+ problem_.temperature(globalPosFace, *eIt));
}
// determine fluid properties at the boundary
Scalar Xw1Bound = BCfluidState.massFrac(wPhaseIdx, wCompIdx);
Scalar Xn1Bound = BCfluidState.massFrac(nPhaseIdx, wCompIdx);
Scalar densityWBound = BCfluidState.density(wPhaseIdx);
- Scalar densityNWBound = BCfluidState.density(nPhaseIdx);
- Scalar viscosityWBound = FluidSystem::phaseViscosity(wPhaseIdx,
- problem_.temperature(globalPosFace, *eIt),
- pressBound, BCfluidState);
- Scalar viscosityNWBound = FluidSystem::phaseViscosity(nPhaseIdx,
- problem_.temperature(globalPosFace, *eIt),
- pressBound+pcBound, BCfluidState);
+ Scalar densityNWBound = BCfluidState.density(nPhaseIdx);
+ Scalar viscosityWBound = FluidSystem::phaseViscosity(wPhaseIdx,
+ problem_.temperature(globalPosFace, *eIt),
+ pressBound, BCfluidState);
+ Scalar viscosityNWBound = FluidSystem::phaseViscosity(nPhaseIdx,
+ problem_.temperature(globalPosFace, *eIt),
+ pressBound+pcBound, BCfluidState);
// average
double densityW_mean = (densityWI + densityWBound) / 2;
@@ -495,9 +495,9 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
K_I.umv(unitDistVec,K);
double velocityW = lambdaW * ((K * unitOuterNormal) * (pressI - pressBound)
- / (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean);
+ / (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean);
double velocityN = lambdaN * ((K * unitOuterNormal) * (pressI - pressBound)
- / (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean);
+ / (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean);
// standardized velocity
double velocityJIw = std::max(-velocityW * faceArea / volume, 0.0);
@@ -534,23 +534,23 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
updFactor[nCompIdx] = - J[1] * faceArea / volume;
// for timestep control
- #define cflIgnoresNeumann
- #ifdef cflIgnoresNeumann
+ #define cflIgnoresNeumann
+ #ifdef cflIgnoresNeumann
factor[0] = 0;
factor[1] = 0;
- #else
+ #else
double inflow = updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW;
if (inflow>0)
- {
- factor[0] = updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW; // =factor in
- factor[1] = -(updFactor[wCompIdx] / densityW /SwmobI + updFactor[nCompIdx] / densityNW / SnmobI); // =factor out
- }
+ {
+ factor[0] = updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW; // =factor in
+ factor[1] = -(updFactor[wCompIdx] / densityW /SwmobI + updFactor[nCompIdx] / densityNW / SnmobI); // =factor out
+ }
else
{
- factor[0] = -(updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW); // =factor in
- factor[1] = updFactor[wCompIdx] / densityW /SwmobI + updFactor[nCompIdx] / densityNW / SnmobI; // =factor out
+ factor[0] = -(updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW); // =factor in
+ factor[1] = updFactor[wCompIdx] / densityW /SwmobI + updFactor[nCompIdx] / densityNW / SnmobI; // =factor out
}
- #endif
+ #endif
}//end neumann boundary
}//end boundary
// correct update Factor by volume error
@@ -569,7 +569,7 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
}// end all intersections
/************************************
- * Handle source term
+ * Handle source term
***********************************/
Dune::FieldVector<double,2> q = problem_.source(globalPos, *eIt);
updateVec[wCompIdx][globalIdxI] += q[wCompIdx];
@@ -577,7 +577,7 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
// account for porosity
sumfactorin = std::max(sumfactorin,sumfactorout)
- / problem_.spatialParameters().porosity(globalPos, *eIt);
+ / problem_.spatialParameters().porosity(globalPos, *eIt);
if ( 1./sumfactorin < dt)
{
diff --git a/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh b/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh
index 1778de4457..10a9c5cceb 100644
--- a/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh
+++ b/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh
@@ -351,7 +351,7 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
}
/******************************************
- * Boundary Face
+ * Boundary Face
******************************************/
if (isIt->boundary())
{
@@ -513,13 +513,13 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
double inflow = updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW;
if (inflow>0)
{
- factor[0] = updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW; // =factor in
- factor[1] = -(updFactor[wCompIdx] / densityW /SwmobI + updFactor[nCompIdx] / densityNW / SnmobI); // =factor out
+ factor[0] = updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW; // =factor in
+ factor[1] = -(updFactor[wCompIdx] / densityW /SwmobI + updFactor[nCompIdx] / densityNW / SnmobI); // =factor out
}
else
{
- factor[0] = -(updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW); // =factor in
- factor[1] = updFactor[wCompIdx] / densityW /SwmobI + updFactor[nCompIdx] / densityNW / SnmobI; // =factor out
+ factor[0] = -(updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW); // =factor in
+ factor[1] = updFactor[wCompIdx] / densityW /SwmobI + updFactor[nCompIdx] / densityNW / SnmobI; // =factor out
}
#endif
}//end neumann boundary
@@ -535,7 +535,7 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
}// end all intersections
/************************************
- * Handle source term
+ * Handle source term
***********************************/
Dune::FieldVector<double,2> q = problem_.source(globalPos, *eIt);
updateVec[wCompIdx][globalIdxI] += q[wCompIdx];
diff --git a/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh b/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh
index 0586824762..3acafbe263 100644
--- a/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh
+++ b/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh
@@ -48,8 +48,8 @@ class PseudoOnePTwoCFluidState : public FluidState<typename GET_PROP_TYPE(TypeTa
public:
- enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
- numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),
+ enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
+ numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),
wPhaseIdx = 0,
nPhaseIdx = 1,};
@@ -135,9 +135,9 @@ public:
// as the mass fractions are calculated, these are used to determine the mole fractions
if (compIdx == wPhaseIdx) //only interested in wComp => X1
{
- moleFrac_ = ( massfracX1_[phaseIdx] / FluidSystem::molarMass(compIdx) ); // = moles of compIdx
+ moleFrac_ = ( massfracX1_[phaseIdx] / FluidSystem::molarMass(compIdx) ); // = moles of compIdx
moleFrac_ /= ( massfracX1_[phaseIdx] / FluidSystem::molarMass(0)
- + (1.-massfracX1_[phaseIdx]) / FluidSystem::molarMass(1) ); // /= total moles in phase
+ + (1.-massfracX1_[phaseIdx]) / FluidSystem::molarMass(1) ); // /= total moles in phase
}
else // interested in nComp => 1-X1
{
@@ -145,7 +145,7 @@ public:
moleFrac_ /= ((1.- massfracX1_[phaseIdx] )/ FluidSystem::molarMass(0)
+ massfracX1_[phaseIdx] / FluidSystem::molarMass(1) ); // /= total moles in phase
}
- return moleFrac_;
+ return moleFrac_;
}
//@}
diff --git a/dumux/decoupled/2p2c/variableclass2p2c.hh b/dumux/decoupled/2p2c/variableclass2p2c.hh
index b1540ed9f5..98dda3ea15 100644
--- a/dumux/decoupled/2p2c/variableclass2p2c.hh
+++ b/dumux/decoupled/2p2c/variableclass2p2c.hh
@@ -88,7 +88,7 @@ private:
ScalarSolutionType saturation_;
PhasePropertyType mobility_; //store lambda for efficiency reasons
ScalarSolutionType capillaryPressure_;
- DimVecElemFaceType velocitySecondPhase_; //necessary??
+ DimVecElemFaceType velocitySecondPhase_; //necessary??
FluidPropertyType density_;
FluidPropertyType numericalDensity_;
@@ -129,7 +129,7 @@ public:
Scalar, dim>& initialVel = *(new Dune::FieldVector<Scalar, dim>(0))) :
VariableClass<TypeTag> (gridView, codim, initialVel), codim_(codim)
{
- initialize2p2cVariables(initialVel, initialSat);
+ initialize2p2cVariables(initialVel, initialSat);
}
//! serialization methods -- same as Dumux::VariableClass2P
@@ -147,7 +147,7 @@ public:
void serializeEntity(std::ostream &outstream, const Element &element)
{
- Dune::dwarn << "here, rather total concentrations should be used!" << std::endl;
+ Dune::dwarn << "here, rather total concentrations should be used!" << std::endl;
int globalIdx = this->elementMapper().map(element);
outstream << this->pressure()[globalIdx] << " " << saturation_[globalIdx];
}
@@ -170,11 +170,11 @@ private:
void initialize2p2cVariables(Dune::FieldVector<Scalar, dim>& initialVel, int initialSat)
{
//resize to grid size
- int size_ = this->gridSize();
- // a) global variables
+ int size_ = this->gridSize();
+ // a) global variables
velocitySecondPhase_.resize(size_);//depends on pressure
velocitySecondPhase_ = initialVel;
- for (int i=0; i<2; i++) //for both phases
+ for (int i=0; i<2; i++) //for both phases
{
density_[i].resize(size_);//depends on pressure
numericalDensity_[i].resize(size_);//depends on pressure
@@ -210,10 +210,10 @@ public:
template<class MultiWriter>
void addOutputVtkFields(MultiWriter &writer)
{
- int size_ = this->gridSize();
+ int size_ = this->gridSize();
if (codim_ == 0)
{
- // pressure & saturation
+ // pressure & saturation
ScalarSolutionType *pressure = writer.template createField<Scalar, 1> (size_);
ScalarSolutionType *saturation = writer.template createField<Scalar, 1> (size_);
*pressure = this->pressure();
diff --git a/dumux/decoupled/common/decoupledproperties.hh b/dumux/decoupled/common/decoupledproperties.hh
index 293d0e5c01..f2f7a327ce 100644
--- a/dumux/decoupled/common/decoupledproperties.hh
+++ b/dumux/decoupled/common/decoupledproperties.hh
@@ -165,10 +165,10 @@ public:
*/
SET_PROP_DEFAULT(TransportSolutionType)
{
- private:
+ private:
typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionType;
- public:
+ public:
typedef typename SolutionType::ScalarSolution type;//!<type for vector of scalar properties
};
diff --git a/dumux/decoupled/common/impet.hh b/dumux/decoupled/common/impet.hh
index ec486f0351..f850ca1c4f 100644
--- a/dumux/decoupled/common/impet.hh
+++ b/dumux/decoupled/common/impet.hh
@@ -101,7 +101,7 @@ public:
*/
void update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec)
{
- // the method is valid for any transported quantity.
+ // the method is valid for any transported quantity.
int transSize = problem.variables().gridSize();
TransportSolutionType transportedQuantity(problem.variables().transportedQuantity());
TransportSolutionType transValueOldIter(problem.variables().transportedQuantity());
diff --git a/test/decoupled/1p/benchmarkresult.hh b/test/decoupled/1p/benchmarkresult.hh
index df5c2b8d4a..ccc98174e7 100644
--- a/test/decoupled/1p/benchmarkresult.hh
+++ b/test/decoupled/1p/benchmarkresult.hh
@@ -36,23 +36,23 @@ namespace Dumux
struct BenchmarkResult
{
private:
- template<int dim>
- struct ElementLayout
- {
- bool contains (Dune::GeometryType gt)
- {
- return gt.dim() == dim;
- }
- };
-
- template<int dim>
- struct FaceLayout
- {
- bool contains (Dune::GeometryType gt)
- {
- return gt.dim() == dim-1;
- }
- };
+ template<int dim>
+ struct ElementLayout
+ {
+ bool contains (Dune::GeometryType gt)
+ {
+ return gt.dim() == dim;
+ }
+ };
+
+ template<int dim>
+ struct FaceLayout
+ {
+ bool contains (Dune::GeometryType gt)
+ {
+ return gt.dim() == dim-1;
+ }
+ };
public:
double relativeL2Error;
@@ -325,23 +325,23 @@ public:
struct ResultEvaluation
{
private:
- template<int dim>
- struct ElementLayout
- {
- bool contains (Dune::GeometryType gt)
- {
- return gt.dim() == dim;
- }
- };
-
- template<int dim>
- struct FaceLayout
- {
- bool contains (Dune::GeometryType gt)
- {
- return gt.dim() == dim-1;
- }
- };
+ template<int dim>
+ struct ElementLayout
+ {
+ bool contains (Dune::GeometryType gt)
+ {
+ return gt.dim() == dim;
+ }
+ };
+
+ template<int dim>
+ struct FaceLayout
+ {
+ bool contains (Dune::GeometryType gt)
+ {
+ return gt.dim() == dim-1;
+ }
+ };
public:
double relativeL2Error;
diff --git a/test/decoupled/2p/test_transport_problem.hh b/test/decoupled/2p/test_transport_problem.hh
index 135b196826..c950ff26f1 100644
--- a/test/decoupled/2p/test_transport_problem.hh
+++ b/test/decoupled/2p/test_transport_problem.hh
@@ -60,49 +60,49 @@ NEW_TYPE_TAG(TransportTestProblem, INHERITS_FROM(DecoupledModel, Transport));
// Set the grid type
SET_PROP(TransportTestProblem, Grid)
{
- typedef Dune::YaspGrid<2> type;
+ typedef Dune::YaspGrid<2> type;
};
// Set the problem property
SET_PROP(TransportTestProblem, Problem)
{
public:
- typedef Dumux::TestTransportProblem<TTAG(TransportTestProblem)> type;
+ typedef Dumux::TestTransportProblem<TTAG(TransportTestProblem)> type;
};
// Set the model properties
SET_PROP(TransportTestProblem, Model)
{
- typedef Dumux::FVSaturation2P<TTAG(TransportTestProblem)> type;
+ typedef Dumux::FVSaturation2P<TTAG(TransportTestProblem)> type;
};
// Set the wetting phase
SET_PROP(TransportTestProblem, WettingPhase)
{
private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
public:
- typedef Dumux::LiquidPhase<Scalar, Dumux::Unit<Scalar> > type;
+ typedef Dumux::LiquidPhase<Scalar, Dumux::Unit<Scalar> > type;
};
// Set the non-wetting phase
SET_PROP(TransportTestProblem, NonwettingPhase)
{
private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
public:
- typedef Dumux::LiquidPhase<Scalar, Dumux::Unit<Scalar> > type;
+ typedef Dumux::LiquidPhase<Scalar, Dumux::Unit<Scalar> > type;
};
// Set the spatial parameters
SET_PROP(TransportTestProblem, SpatialParameters)
{
private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
public:
- typedef Dumux::TestTransportSpatialParams<TypeTag> type;
+ typedef Dumux::TestTransportSpatialParams<TypeTag> type;
};
// Disable gravity
@@ -130,118 +130,118 @@ SET_SCALAR_PROP(TransportTestProblem, CFLFactor, 1.0);
template<class TypeTag = TTAG(TransportTestProblem)>
class TestTransportProblem: public TransportProblem2P<TypeTag, TestTransportProblem<TypeTag> >
{
- typedef TestTransportProblem<TypeTag> ThisType;
- typedef TransportProblem2P<TypeTag, ThisType> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef TestTransportProblem<TypeTag> ThisType;
+ typedef TransportProblem2P<TypeTag, ThisType> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
- enum
- {
- dim = GridView::dimension, dimWorld = GridView::dimensionworld
- };
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
- enum
- {
- wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
- };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
+ };
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GridView::Traits::template Codim<0>::Entity Element;
- typedef typename GridView::Intersection Intersection;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename GridView::Intersection Intersection;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
public:
- TestTransportProblem(const GridView &gridView, const GlobalPosition lowerLeft = 0, const GlobalPosition upperRight = 0) :
- ParentType(gridView), lowerLeft_(lowerLeft), upperRight_(upperRight)
- {
- GlobalPosition vel(0);
- vel[0] = 1e-5;
- this->variables().velocity() = vel;
- this->variables().initializePotentials(vel);
- }
-
- /*!
- * \name Problem parameters
- */
- // \{
-
- /*!
- * \brief The problem name.
- *
- * This is used as a prefix for files generated by the simulation.
- */
- const char *name() const
- {
- return "test_transport";
- }
-
- bool shouldWriteRestartFile() const
- {
- return false;
- }
-
- /*!
- * \brief Returns the temperature within the domain.
- *
- * This problem assumes a temperature of 10 degrees Celsius.
- */
- Scalar temperature(const GlobalPosition& globalPos, const Element& element) const
- {
- return 273.15 + 10; // -> 10°C
- }
-
- // \}
-
-
- Scalar referencePressure(const GlobalPosition& globalPos, const Element& element) const
- {
- return 1e5; // -> 10°C
- }
-
- std::vector<Scalar> source(const GlobalPosition& globalPos, const Element& element)
+ TestTransportProblem(const GridView &gridView, const GlobalPosition lowerLeft = 0, const GlobalPosition upperRight = 0) :
+ ParentType(gridView), lowerLeft_(lowerLeft), upperRight_(upperRight)
+ {
+ GlobalPosition vel(0);
+ vel[0] = 1e-5;
+ this->variables().velocity() = vel;
+ this->variables().initializePotentials(vel);
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const char *name() const
+ {
+ return "test_transport";
+ }
+
+ bool shouldWriteRestartFile() const
+ {
+ return false;
+ }
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature(const GlobalPosition& globalPos, const Element& element) const
+ {
+ return 273.15 + 10; // -> 10°C
+ }
+
+ // \}
+
+
+ Scalar referencePressure(const GlobalPosition& globalPos, const Element& element) const
+ {
+ return 1e5; // -> 10°C
+ }
+
+ std::vector<Scalar> source(const GlobalPosition& globalPos, const Element& element)
{
- return std::vector<Scalar>(2, 0.0);
+ return std::vector<Scalar>(2, 0.0);
}
- BoundaryConditions::Flags bctypeSat(const GlobalPosition& globalPos, const Intersection& intersection) const
- {
- if (globalPos[0] < eps_)
- return Dumux::BoundaryConditions::dirichlet;
- else if (globalPos[0] > upperRight_[0] - eps_)
- return Dumux::BoundaryConditions::outflow;
- else
- return Dumux::BoundaryConditions::neumann;
- }
-
- Scalar dirichletSat(const GlobalPosition& globalPos, const Intersection& intersection) const
- {
- if (globalPos[0] < eps_)
- return 1.0;
- // all other boundaries
- return 0.0;
- }
-
- std::vector<Scalar> neumann(const GlobalPosition& globalPos, const Intersection& intersection) const
- {
- return std::vector<Scalar>(2, 0.0);
- }
-
- Scalar initSat(const GlobalPosition& globalPos, const Element& element) const
- {
- return 0.0;
- }
+ BoundaryConditions::Flags bctypeSat(const GlobalPosition& globalPos, const Intersection& intersection) const
+ {
+ if (globalPos[0] < eps_)
+ return Dumux::BoundaryConditions::dirichlet;
+ else if (globalPos[0] > upperRight_[0] - eps_)
+ return Dumux::BoundaryConditions::outflow;
+ else
+ return Dumux::BoundaryConditions::neumann;
+ }
+
+ Scalar dirichletSat(const GlobalPosition& globalPos, const Intersection& intersection) const
+ {
+ if (globalPos[0] < eps_)
+ return 1.0;
+ // all other boundaries
+ return 0.0;
+ }
+
+ std::vector<Scalar> neumann(const GlobalPosition& globalPos, const Intersection& intersection) const
+ {
+ return std::vector<Scalar>(2, 0.0);
+ }
+
+ Scalar initSat(const GlobalPosition& globalPos, const Element& element) const
+ {
+ return 0.0;
+ }
private:
- GlobalPosition lowerLeft_;
- GlobalPosition upperRight_;
+ GlobalPosition lowerLeft_;
+ GlobalPosition upperRight_;
- static const Scalar eps_ = 1e-6;
+ static const Scalar eps_ = 1e-6;
};
} //end namespace
diff --git a/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh b/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh
index 7289c4eb74..b9dcddf7cd 100644
--- a/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh
+++ b/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh
@@ -67,7 +67,7 @@ public:
const FieldMatrix& intrinsicPermeability (const GlobalPosition& globalPos, const Element& element) const
{
- return constPermeability_;
+ return constPermeability_;
}
double porosity(const GlobalPosition& globalPos, const Element& element) const
diff --git a/test/decoupled/2p2c/test_dec2p2cproblem.hh b/test/decoupled/2p2c/test_dec2p2cproblem.hh
index a45c5f1120..6dcb2f3071 100644
--- a/test/decoupled/2p2c/test_dec2p2cproblem.hh
+++ b/test/decoupled/2p2c/test_dec2p2cproblem.hh
@@ -235,7 +235,7 @@ typename BoundaryConditions::Flags bcTypeTransport(const GlobalPosition& globalP
*/
BoundaryConditions2p2c::Flags bcFormulation(const GlobalPosition& globalPos, const Intersection& intersection) const
{
- return BoundaryConditions2p2c::concentration;
+ return BoundaryConditions2p2c::concentration;
}
//! Value for dirichlet pressure boundary condition \f$ [Pa] \f$.
/*! In case of a dirichlet BC for the pressure equation, the pressure
@@ -266,7 +266,7 @@ Scalar dirichletTransport(const GlobalPosition& globalPos, const Intersection& i
*/
Dune::FieldVector<Scalar,2> neumann(const GlobalPosition& globalPos, const Intersection& intersection) const
{
- Dune::FieldVector<Scalar,2> neumannFlux(0.0);
+ Dune::FieldVector<Scalar,2> neumannFlux(0.0);
// if (globalPos[1] < 15 && globalPos[1]> 5)
// {
// neumannFlux[nPhaseIdx] = -0.015;
diff --git a/test/decoupled/2p2c/test_multiphysics2p2cproblem.hh b/test/decoupled/2p2c/test_multiphysics2p2cproblem.hh
index 607d276869..3d3ac0a6ec 100644
--- a/test/decoupled/2p2c/test_multiphysics2p2cproblem.hh
+++ b/test/decoupled/2p2c/test_multiphysics2p2cproblem.hh
@@ -236,7 +236,7 @@ typename BoundaryConditions::Flags bcTypeTransport(const GlobalPosition& globalP
*/
BoundaryConditions2p2c::Flags bcFormulation(const GlobalPosition& globalPos, const Intersection& intersection) const
{
- return BoundaryConditions2p2c::concentration;
+ return BoundaryConditions2p2c::concentration;
}
/*!
* \copydoc Dumux::TestDecTwoPTwoCProblem::dirichletPress()
@@ -264,7 +264,7 @@ Scalar dirichletTransport(const GlobalPosition& globalPos, const Intersection& i
*/
Dune::FieldVector<Scalar,2> neumann(const GlobalPosition& globalPos, const Intersection& intersection) const
{
- Dune::FieldVector<Scalar,2> neumannFlux(0.0);
+ Dune::FieldVector<Scalar,2> neumannFlux(0.0);
// if (globalPos[1] < 15 && globalPos[1]> 5)
// {
// neumannFlux[nPhaseIdx] = -0.015;
--
GitLab