Hey Shawn,
here, it helps to change the linear solver (or rather the
preconditioner from ILU to SSOR). This unfortunately is a
compile-time choice done in the problem file. I attach my
corresponding problem file, changes are in lines 83 and
86.
Kind regards
Bernd
On Mon, 26 Mar 2012 09:00:18 -0700 (PDT)
Shawn Zhang <[email protected]> wrote:
Hi Bernd,
I switched the permeability between lens and outside,
which gives the non-convergence exception.
I.e.,
permeability = 1e-12
permeabilityLens = 1e-10
Instead of
From: Bernd Flemisch <[email protected]>
To: DuMuX User Forum <[email protected]>
Sent: Monday, March 26, 2012 11:08 AM
Subject: Re: [DuMuX] Permeability heterogeneity and grid
resolution
Hey Shawn,
ok, this is difficult since I cannot reproduce it. Can
you please send the output of your run with 100x100 and
the attached input file? I increased the verbosity of the
linear solver, so maybe we can see more from there.
Kind regards
Bernd
On Mon, 26 Mar 2012 07:53:03 -0700 (PDT)
Shawn Zhang <[email protected]> wrote:
Hi Bernd,
Thanks a lot for your quick attention. However, with
your input file, and grid resolution 100x100, I got the
same exceptions. I also tried to increase
ResidualReduction parameter even further without any
luck.
It puzzles me that the code runs at lower resolution but
not at higher...
-Shawn
From: Bernd Flemisch <[email protected]>
To: DuMuX User Forum <[email protected]>
Sent: Saturday, March 24, 2012 6:07 AM
Subject: Re: [DuMuX] Permeability heterogeneity and grid
resolution
Hey Shawn,
thank you for reporting this. It helps to relax the
threshold for the reduction of the residual that the
linear solver is required to achieve. I attach a
corresponding input file, where the parameter is adjusted
in line 31.
To me, the default threshold of 1e-12 seems to be very
small. We will discuss whether the default value should
be adjusted.
Kind regards
Bernd
On Fri, 23 Mar 2012 12:36:05 -0700 (PDT)
Shawn Zhang <[email protected]> wrote:
Dear All,
I ran into an issue of 1p model not converging. After a
series of tests, I trace it to this canonical situation,
For standard test/boxmodels/1p problem, if I increase
the 2D grid resolution from 10x10 to 100x100, the solver
will not converge, with the following error reported.
This error, however, will go away if I reduce the perm
difference between lens and exterior from 2 order of
magnitude into 1 order of magnitude. This makes me
speculate that the error is related to sharp gradient at
lens/exterior interface caused by the heterogeneity at
high resolution.
I think the reproducer on standard test makes it likely
a serious problem, though I will be very happily
corrected.
Truly appreciate your quick attention.
Best,
-Shawn
root@szdebian:/home/szhang/dumux/dumux-2.1.0/test/boxmodels/1p#
./test_1p -parameterFile test_1p.input
Welcome aboard DuMuX airlines. Please fasten your
seatbelts! Emergency exits are near the time integration.
Initializing problem '1ptest'
Writing result file for "1ptest"
Newton iteration 1 done, relative error = 0.666666
Newton iteration 2 done, relative error = 5.45083e-06
Solve: M deltax^k = rNewton: Caught exception:
"NumericalProblem
[newtonSolveLinear:../../../dumux/nonlinear/newtoncontroller.hh:398]:
Linear solver did not converge"
Newton solver did not converge with dt=1 seconds.
Retrying with time step of 0.5 seconds
Newton iteration 1 done, relative error = 0.666666
Newton iteration 2 done, relative error = 5.45083e-06
Solve: M deltax^k = rNewton: Caught exception:
"NumericalProblem
[newtonSolveLinear:../../../dumux/nonlinear/newtoncontroller.hh:398]:
Linear solver did not converge"
....
Newton solver did not converge with dt=0.0078125
seconds. Retrying with time step of 0.00390625 seconds
Newton iteration 1 done, relative error = 0.666666
Newton iteration 2 done, relative error = 5.45083e-06
Solve: M deltax^k = rNewton: Caught exception:
"NumericalProblem
[newtonSolveLinear:../../../dumux/nonlinear/newtoncontroller.hh:398]:
Linear solver did not converge"
Newton solver did not converge with dt=0.00390625
seconds. Retrying with time step of 0.00195312 seconds
Newton iteration 1 done, relative error = 0.666666
Newton iteration 2 done, relative error = 5.45083e-06
Solve: M deltax^k = rNewton: Caught exception:
"NumericalProblem
[newtonSolveLinear:../../../dumux/nonlinear/newtoncontroller.hh:398]:
Linear solver did not converge"
Newton solver did not converge with dt=0.00195312
seconds. Retrying with time step of 0.000976562 seconds
Dune reported error: Dune::MathError
[timeIntegration:../../../dumux/boxmodels/common/boxproblem.hh:485]:
Newton solver didn't converge after 10 time-step
divisions. dt=0.000976562
___________________________________________________________
Bernd Flemisch phone: +49 711 685
69162
IWS, Universitaet Stuttgart fax: +49 711 685
67020
Pfaffenwaldring 61 email:
[email protected]
D-70569 Stuttgart url:
www.hydrosys.uni-stuttgart.de
___________________________________________________________
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
___________________________________________________________
Bernd Flemisch phone: +49 711 685
69162
IWS, Universitaet Stuttgart fax: +49 711 685
67020
Pfaffenwaldring 61 email:
[email protected]
D-70569 Stuttgart url:
www.hydrosys.uni-stuttgart.de
___________________________________________________________
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
permeability = 1e-10
permeabilityLens = 1e-12
Sorry about the ignorance. I should notice this the
first time. Can you give a try and see whether you can
reproduce?
-Shawn
___________________________________________________________
Bernd Flemisch phone: +49 711 685
69162
IWS, Universitaet Stuttgart fax: +49 711 685
67020
Pfaffenwaldring 61 email:
[email protected]
D-70569 Stuttgart url:
www.hydrosys.uni-stuttgart.de
___________________________________________________________
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2009 by Onur Dogan *
* Copyright (C) 2009 by Andreas Lauser *
* Institute for Modelling Hydraulic and Environmental Systems *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief A test problem for the one-phase box model:
* water is flowing from bottom to top through and around a low permeable lens.
*/
#ifndef DUMUX_1PTEST_PROBLEM_HH
#define DUMUX_1PTEST_PROBLEM_HH
#if HAVE_UG
#include <dune/grid/io/file/dgfparser/dgfug.hh>
#endif
#if HAVE_ALUGRID
#include <dune/grid/io/file/dgfparser/dgfalu.hh>
#endif
#include <dune/grid/io/file/dgfparser/dgfs.hh>
#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
#include <dune/grid/io/file/gmshreader.hh>
#include <dumux/boxmodels/1p/1pmodel.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/fluidsystems/liquidphase.hh>
#include "1ptestspatialparameters.hh"
namespace Dumux
{
template <class TypeTag>
class OnePTestProblem;
namespace Properties
{
NEW_TYPE_TAG(OnePTestProblem, INHERITS_FROM(BoxOneP));
SET_PROP(OnePTestProblem, Fluid)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleH2O<Scalar> > type;
};
// Set the grid type
SET_PROP(OnePTestProblem, Grid)
{
typedef Dune::SGrid<2, 2> type;
//typedef Dune::YaspGrid<2> type;
//typedef Dune::UGGrid<2> type;
//typedef Dune::ALUSimplexGrid<2,2> type;
};
// Set the problem property
SET_PROP(OnePTestProblem, Problem)
{ typedef Dumux::OnePTestProblem<TypeTag> type; };
// Set the spatial parameters
SET_PROP(OnePTestProblem, SpatialParameters)
{ typedef Dumux::OnePTestSpatialParameters<TypeTag> type; };
// Linear solver settings
SET_TYPE_PROP(OnePTestProblem, LinearSolver, Dumux::BoxCGSSORSolver<TypeTag> );
SET_INT_PROP(OnePTestProblem, LinearSolverVerbosity, 0);
SET_SCALAR_PROP(OnePTestProblem, LinearSolverResidualReduction, 1e-12);
SET_INT_PROP(OnePTestProblem, PreconditionerIterations, 3);
SET_SCALAR_PROP(OnePTestProblem, PreconditionerRelaxation, 1.0);
// Enable gravity
SET_BOOL_PROP(OnePTestProblem, EnableGravity, true);
}
/*!
* \ingroup OnePBoxModel
* \ingroup BoxTestProblems
* \brief Test problem for the one-phase box model:
* water is flowing from bottom to top through and around a low permeable lens.
*
* The domain is box shaped. All sides are closed (Neumann 0 boundary)
* except the top and bottom boundaries (Dirichlet), where water is
* flowing from bottom to top.
*
* In the middle of the domain, a lens with low permeability (\f$K=10e-12\f$)
* compared to the surrounding material (\f$ K=10e-10\f$) is defined.
*
* To run the simulation execute the following line in shell:
* <tt>./test_1p -parameterFile test_1p.input</tt>
* The same parameter file can be also used for 3d simulation but you need to change line
* <tt>typedef Dune::SGrid<2,2> type;</tt> to
* <tt>typedef Dune::SGrid<3,3> type;</tt> in the problem file
* and use <tt>1p_3d.dgf</tt> in the parameter file.
*/
template <class TypeTag>
class OnePTestProblem : public OnePBoxProblem<TypeTag>
{
typedef OnePBoxProblem<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
// copy some indices for convenience
typedef typename GET_PROP_TYPE(TypeTag, OnePIndices) Indices;
enum {
// Grid and world dimension
dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
// indices of the primary variables
pressureIdx = Indices::pressureIdx
};
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, 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, FVElementGeometry) FVElementGeometry;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
public:
OnePTestProblem(TimeManager &timeManager, const GridView &gridView)
: ParentType(timeManager, gridView)
{
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \brief The problem name.
*
* This is used as a prefix for files generated by the simulation.
*/
const char *name() const
{ return "1ptest"; }
/*!
* \brief Return the temperature within the domain.
*
* This problem assumes a temperature of 10 degrees Celsius.
*/
Scalar temperature() const
{ return 273.15 + 10; } // 10C
void sourceAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
values = 0;
}
// \}
/*!
* \name Boundary conditions
*/
// \{
/*!
* \brief Specify which kind of boundary condition should be
* used for which equation on a given boundary segment.
*/
void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
{
const GlobalPosition globalPos = vertex.geometry().center();
double eps = 1.0e-3;
if (globalPos[dim-1] < eps || globalPos[dim-1] > this->bboxMax()[dim-1] - eps)
values.setAllDirichlet();
else
values.setAllNeumann();
}
/*!
* \brief Evaluate the boundary conditions for a Dirichlet
* boundary segment.
*
* For this method, the \a values parameter stores primary variables.
*/
void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
{
double eps = 1.0e-3;
const GlobalPosition globalPos = vertex.geometry().center();
if (globalPos[dim-1] < eps) {
values[pressureIdx] = 2.0e+5;
}
else if (globalPos[dim-1] > this->bboxMax()[dim-1] - eps) {
values[pressureIdx] = 1.0e+5;
}
}
/*!
* \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 component. Negative values mean
* influx.
*/
using ParentType::neumann;
void neumann(PrimaryVariables &values,
const Element &element,
const FVElementGeometry &fvElemGeom,
const Intersection &is,
int scvIdx,
int boundaryFaceIdx) const
{
// const GlobalPosition &globalPos = fvElemGeom.boundaryFace[boundaryFaceIdx].ipGlobal;
values[pressureIdx] = 0;
}
// \}
/*!
* \name Volume terms
*/
// \{
/*!
* \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);
values[pressureIdx] = 1.0e+5;// + 9.81*1.23*(20-globalPos[dim-1]);
}
// \}
};
} //end namespace
#endif
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
