Hey Shawn, On Mon, 2 Apr 2012 06:35:49 -0700 (PDT) Shawn Zhang <[email protected]> wrote:
No. Pretty much anything can happen. Most often the uninitialized z-values will be zero, thus there will be no lens. But also other numbers can be in there, making it unpredictable. So you really have to incorporate the corresponding stuff from my 1ptestspatialparameters.hh file.Hi Bernd, Thanks for the suggestions.I am using the release. All 3D runs work if the resolution is low (say, 10x10x10), but diverge for higher resolution. If the problem is related to lens z-coordinate being undefined, the code should break regardless the resolution, do you agree?
Moreover, in the example I attached in my last email, I assigned a tensor permeability for each cell. Can you please take a quick look, and suggest where I did wrong?
This seems fine.
Again, the code runs with lower resolution, but not higher.No need to become nervous. I attach my current problem file and input file with which it works in 3d for me up to 50x50x50, I did not test larger.Waiting nervously for your advice.
This time, the trick for me was to use BiCGStabSSOR instead of CGSSOR as a solver. And I even can come up with a handwaving explanation. In principle, the matrix should be symmetric. But this is not explicitly enforced in the code, so a_ij and a_ji are calculated separately. Due to numerical errors, it will be slightly unsymmetric and the CG solver gets into trouble. BiCGStab does not care.
Hope this helps. Kind regards Bernd
Best, -Shawn From: Bernd Flemisch <[email protected]>To: DuMuX User Forum <[email protected]> Sent: Monday, April 2, 2012 3:17 AM Subject: Re: [DuMuX] Permeability heterogeneity and grid resolutionHey Shawn,thank you for your further investigations. You are right that the options for the linear solver are not well documented yet. You can choose one out of dumux/linear/boxlinearsolver.hh. The parameters are documented in dumux/linear/linearsolverproperties.hh. About the actual meaning of, e.g., a relaxation parameter or a preconditioner, you should consider your favorite linear algebra textbook.Concerning your 3d calculations: are you using the release or the trunk? If you use the release, 3d can break, since the z-coordinate of the lens is not set. You can use the attached 1ptestspatialparameters.hh. This is already fixed in the trunk and will be fixed in a bugfix release 2.1.1.Kind regards Bernd On Sat, 31 Mar 2012 19:02:42 -0700 (PDT) Shawn Zhang <[email protected]> wrote:Hi Bernd and All,I confirm that Bernd's input file resolved the convergence problem I had. Thanks.However (I hate it when people do this to me, but sorry I can't find other way out), I am still puzzled by this convergence issue. It seems the following factors are influencing the convergence,1. Resolution. The higher the resolution, the more difficult the convergence. 2. Heterogeneity. The higher the heterogeneity, the more difficult the convergence. 3. Linear solver type. Not sure what are the options and how they are related to convergence. 4. Solve controls: relaxation, preconditioners, and number of iterations. I only understand the last one.I most puzzled by #1. Why higher grid resolution causes non-convergence? Attached a more specific problem I am working on. The set up converge at 10x10x10, but diverge at 50x50x50 or higher. Would you mind taking a look and suggest how to make it converge? I need to run it at 200x200x200 resolution.In addition to the problem, I also modified 1pmodel.hh file to handle tensor perm. I also have a question about why using Scalar tensor in all the boxmodels tests, but this is not critical (unless it is related to the convergence issue).Sorry about all the questions, but I am about to pull my hair now.-Shawn From: Bernd Flemisch <[email protected]>To: DuMuX User Forum <[email protected]> Sent: Thursday, March 29, 2012 5:23 AM Subject: Re: [DuMuX] Permeability heterogeneity and grid resolutionSorry, wrong file. Attached is the correct one. Kind regards On Thu, 29 Mar 2012 10:36:57 +0200 Bernd Flemisch <[email protected]> wrote:Hey Shawn,I was not too specific about what to change: you should set PreconditionerIterations to 3, not PreconditionerRelaxation. You can also set this in the input file, it overrides the settings of the problem file.I attach my current working input file for your convenience.Kind regards Bernd On Wed, 28 Mar 2012 12:42:13 -0700 (PDT) Shawn Zhang <[email protected]> wrote:Hi Bernd and All,Sorry if I appear to be sticky about the convergence issue, but it kept coming back and bite me.Coming back to this test problem, with the attached 1p project, I have the following issue.Standard test: worked.Increase grid size to 100x100, not converging. It worked after I, increase # iteration to 1000, reduce convergence criterion to 10e-6. Increase grid size to 1000x1000, not converging. I changed the linear solver from Dumux::BoxCGILU0Solver to Dumux::BoxCGSSORSolver, and PreconditionerRelaxation from 1.0 to 3.0 as suggested by Bernd, it is still not converging...I attached the project folder.Do you think this maybe related to the Oracle VBox based Debian installation?Thanks. -Shawn From: Bernd Flemisch <[email protected]>To: DuMuX User Forum <[email protected]> Sent: Tuesday, March 27, 2012 4:01 AM Subject: Re: [DuMuX] Permeability heterogeneity and grid resolutionHey Shawn,with the attached input file, it works for 1000x1000 and four orders difference. I increased the maximum number of iterations to 1000 and the residual reduction to 1e-6.Kind regards Bernd On Mon, 26 Mar 2012 11:44:21 -0700 (PDT) Shawn Zhang <[email protected]> wrote:Well, I celebrated too fast.My goal is to determine, in a system of two material with different perm, when one material can be treated as having negligible impact on pressure drop.The two parameters influencing this are, 1. Permeability difference between two materials. 2. Grid resolution.I use the 1p with lens test problem, and increasing the perm of the lens material.At 100x100 resolution, I derive the attached plot (pressure plot along y=0.5 line), which is saying that when the permeability difference is greater than 3 order of magnitude, the material with high permeability can be seen as open channel flow comparing to the low permeable phase, and the pressure drop is neglegible. Moreover, further increase of the perm difference will have a neglegible effect.At 10x10 resolution, we start seeing this when the perm difference is 2 order of magnitude.When I increase the resolution to 1000x1000, with the permeability difference 4 order of magnitude, the solver does not converge again.Is this the expected behavior, as I am trying to solve large discontinuity caused by the big perm difference?-Shawn From: Bernd Flemisch <[email protected]>To: DuMuX User Forum <[email protected]> Sent: Monday, March 26, 2012 12:45 PM Subject: Re: [DuMuX] Permeability heterogeneity and grid resolutionHey 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 resolutionHey 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 resolutionHey 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, -Shawnroot@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-06Solve: 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 secondsNewton iteration 1 done, relative error = 0.666666 Newton iteration 2 done, relative error = 5.45083e-06Solve: 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 secondsNewton iteration 1 done, relative error = 0.666666 Newton iteration 2 done, relative error = 5.45083e-06Solve: 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 secondsNewton iteration 1 done, relative error = 0.666666 Newton iteration 2 done, relative error = 5.45083e-06Solve: 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-12Sorry 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___________________________________________________________ _______________________________________________ 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___________________________________________________________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______________________________________________________________________________________________________________________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
___________________________________________________________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
___________________________________________________________
test_1p.input
Description: Binary data
// -*- 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<3, 3> 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::BoxBiCGStabSSORSolver<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
