[test][2p][seq] Update 2p lense problem

test/porousmediumflow/2p/sequential/impes/pressureproblem.hh

@@ -59,7 +59,7 @@ SET_PROP(TwoPImpes, FluidSystem)
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
 #if PROBLEM == 2
-    using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Trichloroethene<Scalar> >;
+    using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
 #else
     using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> >;
 #endif
@@ -135,11 +135,11 @@ public:
     BoundaryTypes values;
     if (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos))
     {
-        values.setAllDirichlet();
+        values.setAllNeumann();
     }
     else
     {
-        values.setAllNeumann();
+        values.setAllDirichlet();
     }
     return values;
 #else
@@ -169,12 +169,8 @@ public:
         Scalar pRef = 1.0e5;
 
         using WettingPhase = typename GET_PROP(TypeTag, FluidSystem)::WettingPhase;
-        if (onUpperBoundary_(globalPos))
-        {
-            pRef = 2.0e5;
-        }
-
-        values[pressureIdx] = (pRef - (this->fvGridGeometry().bBoxMax()- globalPos) * this->gravity() * WettingPhase::density(temperature(), pRef));
+        values[pressureIdx] = (pRef - (this->fvGridGeometry().bBoxMax()- globalPos)
+                                       * this->gravity() * WettingPhase::density(temperature(), pRef));
 
         return values;
     }
@@ -192,9 +188,10 @@ public:
      */
     NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
+        using NonwettingPhase = typename GET_PROP(TypeTag, FluidSystem)::NonwettingPhase;
         NumEqVector values(0.0);
-        //if(isInlet(globalPos))
-        //    values[pressureEqIdx] = -inFlux_;
+        if(isInlet(globalPos))
+            values[pressureEqIdx] = -inFlux_/NonwettingPhase::density(temperature(), 1.0e5);
 
         return values;
     }

test/porousmediumflow/2p/sequential/impes/test_impes.input

 [TimeLoop]
 DtInitial = 1e2# [s]
-MaxTimeStepSize = 1
+MaxTimeStepSize = 10
 TEnd = 5e4 # [s]
+OutputTimeInterval = 10 # [s]
 
 [Grid]
 #File = ../grids/test_mpfa2p.dgf #grid for buckley-leverett or mcwhorter problem
-#Refinement = 0
+Refinement = 1
 LowerLeft = 0 0
 UpperRight = 20 10
-Cells = 100 50
+Cells = 10 10
 
 [SpatialParams]
 BackgroundPermeability = 1e-10
@@ -28,7 +29,7 @@ LenseEntryPressure = 0
 
 [Problem]
 Name = test_lenses
-EnableGravity = true
+EnableGravity = false
 
 InletWidth = 2
 InjectionFlux = 0.1
@@ -36,3 +37,9 @@ InjectionFlux = 0.1
 [Component]
 LiquidDensity = 1000
 LiquidKinematicViscosity = 1e-6
+
+[Newton]
+MaxRelativeShift = 1.0e-14
+
+[Assembly]
+NumericDifference.BaseEpsilon =  1.0e-12

test/porousmediumflow/2p/sequential/impes/test_impes_fv.cc

@@ -170,6 +170,13 @@ int main(int argc, char** argv) try
     VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
     vtkWriter.write(0.0);
 
+    using VtkOutputFieldsTransport = typename GET_PROP_TYPE(TwoPTransportTT, VtkOutputFields);
+
+    std::string nameTransport = twoPImpesProblem->name() + "_transport";
+    VtkOutputModule<GVTwoPTransport, SVTwoPTransport> vtkWriterTransport(*gridVarTwoPTransport, x_s, nameTransport, "",Dune::VTK::conforming);
+    VtkOutputFieldsTransport::init(vtkWriterTransport); //!< Add model specific output fields
+    vtkWriterTransport.write(0.0);
+
     // instantiate time loop
     auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
     timeLoop->setMaxTimeStepSize(maxDt);
@@ -201,7 +208,11 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>();
 
     //! set some check points for the time loop
-    timeLoop->setPeriodicCheckPoint(4);
+    timeLoop->setPeriodicCheckPoint(getParam<Scalar>("TimeLoop.OutputTimeInterval"));
+
+    // the non-linear solver
+    using NewtonSolver = Dumux::NewtonSolver<TwoPImpesAssembler, LinearSolver>;
+    NewtonSolver nonLinearSolver(twoPImpesAssembler, linearSolver);
 
     //! start the time loop
     timeLoop->start();
@@ -237,16 +248,19 @@ int main(int argc, char** argv) try
 
         twoPImpesProblem->spatialParams().setWettingSaturation(satVals);
 
-        // set previous solution for storage evaluations
+//        // set previous solution for storage evaluations
         twoPImpesAssembler->setPreviousSolution(x_p_old);
-
-        twoPImpesAssembler->assembleJacobianAndResidual(x_p);
-
-        SVTwoPImpes xpDelta(x_p);
-        linearSolver->solve(*Ap, xpDelta, *rp);
-
-        x_p -= xpDelta;
-        gridVarTwoPImpes->update(x_p);
+//
+//        twoPImpesAssembler->assembleJacobianAndResidual(x_p);
+//
+//        SVTwoPImpes xpDelta(x_p);
+//        linearSolver->solve(*Ap, xpDelta, *rp);
+//
+//        x_p -= xpDelta;
+//        gridVarTwoPImpes->update(x_p);
+
+        // solve the non-linear system with time step control
+        nonLinearSolver.solve(x_p, *timeLoop);
 
         // make the new solution the old solution
         x_p_old = x_p;
@@ -256,7 +270,7 @@ int main(int argc, char** argv) try
         std::vector<Scalar> volumeFlux(twoPImpesFvGridGeometry->numScvf(), 0.0);
 
         using FluxVariables =  typename GET_PROP_TYPE(TwoPImpesTT, FluxVariables);
-        auto upwindTerm = [](const auto& volVars) { return volVars.mobility(0); };
+        auto upwindTerm = [](const auto& volVars) { return 1.0; };
         for (const auto& element : elements(leafGridView))
         {
             auto fvGeometry = localView(*twoPImpesFvGridGeometry);
@@ -320,7 +334,10 @@ int main(int argc, char** argv) try
 
         // write vtk output on check points
         if (timeLoop->isCheckPoint() || timeLoop->finished())
+        {
             vtkWriter.write(timeLoop->time());
+            vtkWriterTransport.write(timeLoop->time());
+        }
 
         // set new dt
         timeLoop->setTimeStepSize(dt);

test/porousmediumflow/2p/sequential/impes/transportproblem.hh

@@ -58,7 +58,7 @@ SET_PROP(TwoPTransport, FluidSystem)
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
 #if PROBLEM == 2
-    using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Trichloroethene<Scalar> >;
+    using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
 #else
     using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> >;
 #endif
@@ -134,11 +134,11 @@ public:
         BoundaryTypes values;
         if (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos))
         {
-            values.setAllDirichlet();
+            values.setAllNeumann();
         }
         else
         {
-            values.setAllNeumann();
+            values.setAllDirichlet();
         }
         return values;
 #else
@@ -162,17 +162,8 @@ public:
     PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
     {
         PrimaryVariables values;
-//        if (onLeftBoundary_(globalPos))
-//                values[saturationIdx] = 0.8;
-
-        using WettingPhase = typename GET_PROP(TypeTag, FluidSystem)::WettingPhase;
         values[saturationIdx] = 1.0;
 
-        if (isInlet(globalPos))
-        {
-            values[saturationIdx] = 0.0;
-        }
-
         return values;
     }
 

test/porousmediumflow/2p/sequential/impes/transportspatialparams.hh

@@ -165,6 +165,9 @@ public:
                       const ElementVolumeVariables& elemVolVars,
                       const SubControlVolumeFace& scvf) const
     {
+//        GlobalPosition vel(0.0);
+//        vel[dimWorld-1] = -1.0e-6;
+//        return vel*scvf.unitOuterNormal();
         return volumeFlux_[scvf.index()];
     }