GEOS-DEV
GMRES/MGR is stalling, UMFPACK is out of memory for >=32k mesh in SPE10 mechanical extension model
Hi everybody! I am comparing thermoporoelastic modeling in GEOS (FEM) and DARTS (FVM) in the case of SPE10 mechanical extension model introduced in Garipov et al, 2018. This is an extension of SPE10 dataset to geomechanical model, where bulk modulus becomes heterogeneous (linear against porosity) along with porosity and permeability. I use some volume-averaging to upscale the model to various resolutions. By now, i am comparing only initialization, where initial pressure and temperature are defined and i expect to see the effect of applied boundary conditions (uniform load from top, rollers at sides and bottom). I build up poroelastic and thermoporoelastic models of different resolutions (10x10x10, 20x40x40, 40x80x40).
The problem i am experiencing is that the run (specifically singlePhasePoroelasticityPreEquilibrationStep) is stalling for thermoporoelastic setup for moderate and high grid resolutions (20x40x40=32k, 40x80x40=128k) while running smoothly for coarse grid (10x10x10) and for all grids in pororoelastic setup. I tried two solvers:
<LinearSolverParameters directParallel="0"/><LinearSolverParameters solverType="gmres" preconditionerType="mgr"/>The first one (direct) run out of memory for 32k grid! See output below. The iterative one is stalling like this:
(base) bash-4.4$ ./../../../GEOS/build-environment-release/bin/geosx -i thermoporoelastic_fim.xml
Num ranks: 1
GEOS version: 0.2.0 (develop, sha1: c11ad8ecd)
- c++ compiler: gcc 11.3.0
- MPI version: Intel(R) MPI Library 2021.9 for Linux* OS
- HDF5 version: 1.12.1
- Conduit version: 0.8.2
- VTK version: 9.2.6
- RAJA version: 2023.6.1
- umpire version: 2023.6.0
- adiak version: ..
- caliper version: 2.10.0
- METIS version: METIS_VERSION
- PARAMETIS version: PARAMETIS_VERSION
- scotch version: 7.0.3
- superlu_dist version: 6.3.0
- suitesparse version: 5.7.9
- hypre version: v2.30.0-44-geab5f9f7f (master)
- Python3 version: 3.10.10
Started at 2024-04-02 12:27:02.476506410
Opened XML file: /oahu/data/avnovikov/geos/my_models/SPE10_mech/geos/thermoporoelastic_fim.xml
Included additionnal XML file: /oahu/data/avnovikov/geos/my_models/SPE10_mech/geos/thermoporoelastic_base.xml
Adding Solver of type SinglePhasePoromechanicsReservoir, named reservoirSystem
Adding Solver of type SinglePhasePoromechanics, named singlePhasePoroelasticity
Adding Solver of type SolidMechanicsLagrangianSSLE, named linearElasticity
Adding Solver of type SinglePhaseFVM, named singlePhaseFlow
Adding Solver of type SinglePhaseWell, named singlePhaseWell
Adding Event: PeriodicEvent, outputs
Adding Event: SoloEvent, singlePhasePoroelasticityPreEquilibrationStep
Adding Event: PeriodicEvent, solverApplicationsEquilibration
Adding Event: SoloEvent, singlePhasePoroelasticityPostEquilibrationStep
Adding Event: PeriodicEvent, linearElasticityStatistics
Adding Event: PeriodicEvent, singlePhaseFlowStatistics
Adding Event: PeriodicEvent, solverApplications
Adding Event: PeriodicEvent, restarts
Adding Output: VTK, vtkOutput
Adding Output: Restart, restartOutput
Adding Mesh: VTKMesh, mesh
Adding Geometric Object: Box, xneg
Adding Geometric Object: Box, xpos
Adding Geometric Object: Box, yneg
Adding Geometric Object: Box, ypos
Adding Geometric Object: Box, zneg
Adding Geometric Object: Box, zpos
TableFunction: temperature0
TableFunction: temperature_side
TableFunction: pressure_side
Adding Object CellElementRegion named Domain from ObjectManager::Catalog.
VTK `vtkOutput (thermoporoelastic_fim.xml, l.130)`: found 1 fields to plot in `fieldNames`, in addition to all fields with `plotLevel` smaller or equal to 1.
reservoirSystem: found SinglePhasePoromechanics solver named singlePhasePoroelasticity
reservoirSystem: found SinglePhaseWell solver named singlePhaseWell
singlePhasePoroelasticity: found SolverBase solver named singlePhaseFlow
singlePhasePoroelasticity: found SolidMechanics_LagrangianFEM solver named linearElasticity
VTKMesh 'mesh': reading mesh from /oahu/data/avnovikov/geos/my_models/SPE10_mech/geos/input/data_20_40_40/geos_input.vtu
reading the dataset...
redistributing mesh...
Generating global Ids from VTK mesh
Rank 0: Local mesh size: 32000
finding neighbor ranks...
done!
VTKMesh 'mesh': generating GEOSX mesh data structure
preprocessing...
writing nodes...
writing cells...
Importing cell block 99991_hexahedra
writing surfaces...
building connectivity maps...
done!
Number of nodes: 35301
Number of elems: 32000
C3D8: 32000
Load balancing: min avg max
(element/rank): 32000 32000 32000
regionQuadrature: meshBodyName, meshLevelName, regionName, subRegionName = mesh, Level0, Domain, 99991_hexahedra
mesh/Level0/Domain/99991_hexahedra/fluid allocated 8 quadrature points
mesh/Level0/Domain/99991_hexahedra/porousRock allocated 8 quadrature points
mesh/Level0/Domain/99991_hexahedra/thermalCond allocated 8 quadrature points
mesh: importing field data from mesh dataset
Importing field biotCoefficient into rockPorosity_biotCoefficient on Domain/99991_hexahedra
Importing field bulkModulus into skeleton_bulkModulus on Domain/99991_hexahedra
Importing field permeability into rockPerm_permeability on Domain/99991_hexahedra
Importing field porosity into rockPorosity_referencePorosity on Domain/99991_hexahedra
Importing field biotCoefficient into rockPorosity_biotCoefficient on Domain/99991_hexahedra
Importing field bulkModulus into skeleton_bulkModulus on Domain/99991_hexahedra
Importing field permeability into rockPerm_permeability on Domain/99991_hexahedra
Importing field porosity into rockPorosity_referencePorosity on Domain/99991_hexahedra
Importing field biotCoefficient into rockPorosity_biotCoefficient on Domain/99991_hexahedra
Importing field bulkModulus into skeleton_bulkModulus on Domain/99991_hexahedra
Importing field permeability into rockPerm_permeability on Domain/99991_hexahedra
Importing field porosity into rockPorosity_referencePorosity on Domain/99991_hexahedra
TableFunction: equil_99991_hexahedra_table
------------------- TIMESTEP START -------------------
- Time: -3168y, -319d, -4h-1m-4s (-100000000000 s)
- Delta Time: 3168y, 319d, 04h01m04s (100000000000 s)
- Cycle: 0
------------------------------------------------------
Time: -1.00e+11 s, dt: 100000000000 s, Cycle: 0
Task `singlePhasePoroelasticityPreEquilibrationStep`: at time -100000000000s, physics solver `singlePhasePoroelasticity` is set to perform stress initialization during the next time step(s)
Task `singlePhasePoroelasticityPreEquilibrationStep`: at time -100000000000s, physics solver `linearElasticity` is resetting total displacement and velocity to zero
Attempt: 0, ConfigurationIter: 0, NewtonIter: 0
( Rflow ) = ( 0.00e+00 ) ( Renergy ) = ( 0.00e+00 ) ( Rsolid ) = ( 1.72e+02 ) ( Rwell ) = ( 0.00e+00 ) ( R ) = ( 1.72e+02 )
I checked input porosities, permeabilities and bulkModulus in 20x40x40 mesh - they remain within reasonable limits:
- (0.0400641, 0.32541) for porosity
- (1983405555.5555558, 7410333333.333334) for bulkModulus
- (1.005522803e-21, 1.2564223900000002e-11) for permeability
- Biot coefficient is uniform = 1.0
Stacktrace with GMRES/MGR after Ctrl+C is following:
** StackTrace of 16 frames **
Frame 0: /lib64/libc.so.6
Frame 1: hypre_GaussElimSetup
Frame 2: hypre_MGRSetup
Frame 3: geos::HyprePreconditioner::setup(geos::HypreMatrix const&)
Frame 4: geos::HypreSolver::setup(geos::HypreMatrix const&)
Frame 5: geos::SolverBase::solveLinearSystem(geos::DofManager const&, geos::HypreMatrix&, geos::HypreVector&, geos::HypreVector&)
Frame 6: geos::SolverBase::solveNonlinearSystem(double const&, double const&, int, geos::DomainPartition&)
Frame 7: geos::SolverBase::nonlinearImplicitStep(double const&, double const&, int, geos::DomainPartition&)
Frame 8: geos::SolverBase::solverStep(double const&, double const&, int, geos::DomainPartition&)
Frame 9: geos::CoupledSolver<geos::SinglePhasePoromechanics<geos::SinglePhaseBase>, geos::SinglePhaseWell>::solverStep(double const&, double const&, int, geos::DomainPartition&)
Frame 10: geos::SolverBase::execute(double, double, int, int, double, geos::DomainPartition&)
Frame 11: geos::EventBase::execute(double, double, int, int, double, geos::DomainPartition&)
Frame 12: geos::EventManager::run(geos::DomainPartition&)
Frame 13: geos::GeosxState::run()
Frame 14: main
Frame 15: __libc_start_main
Frame 16: _start
=====
Stacktrace with direct solver after Ctrl+C is following:
UMFPACK V5.7.9 (Oct 20, 2019): ERROR: out of memory
***** ERROR
***** LOCATION: /oahu/data/avnovikov/geos/GEOS/src/coreComponents/linearAlgebra/interfaces/direct/SuiteSparse.cpp:158
***** Controlling expression (should be false): true
***** Rank 0: SuiteSparse: umfpack_dl_numeric failed.
** StackTrace of 12 frames **
Frame 0: geos::SuiteSparse<geos::HypreInterface>::setup(geos::HypreMatrix const&)
Frame 1: geos::SolverBase::solveLinearSystem(geos::DofManager const&, geos::HypreMatrix&, geos::HypreVector&, geos::HypreVector&)
Frame 2: geos::SolverBase::solveNonlinearSystem(double const&, double const&, int, geos::DomainPartition&)
Frame 3: geos::SolverBase::nonlinearImplicitStep(double const&, double const&, int, geos::DomainPartition&)
Frame 4: geos::SolverBase::solverStep(double const&, double const&, int, geos::DomainPartition&)
Frame 5: geos::CoupledSolver<geos::SinglePhasePoromechanics<geos::SinglePhaseBase>, geos::SinglePhaseWell>::solverStep(double const&, double const&, int, geos::DomainPartition&)
Frame 6: geos::SolverBase::execute(double, double, int, int, double, geos::DomainPartition&)
Frame 7: geos::EventBase::execute(double, double, int, int, double, geos::DomainPartition&)
Frame 8: geos::EventManager::run(geos::DomainPartition&)
Frame 9: geos::GeosxState::run()
Frame 10: main
Frame 11: __libc_start_main
Frame 12: _start
=====
Abort(1) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0
Debugging does not show anything valuable so far.
To Reproduce Could anybody please just clone and try to run it?
git clone https://github.com/av-novikov/spe10_mech_geos.gitcd spe10_mech_geos./path/to/geosx -i thermoporoelastic_fim.xmlto run thermoporoelastic model and see stalled run./path/to/geosx -i poroelastic_fim.xmlto run poroelastic model and see model is running- Adjust
file="input/data_20_40_40/geos_input.vtu"/>tofile="input/data_10_10_10/geos_input.vtu"/>in thermoporoelastic_base.xml to make thermoporoelastic_fim.xml running without being stalled.
Expected behavior Should run smoothly with GMRES/MGR linear solver. No stalling as in poroelastic example.
Platform (please complete the following information):
- oahu.citg.tudelft.nl, RHEL 8.8 (Ootpa), Linux 4.18.0-477.13.1.el8_8.x86_64
- 639Gi available unlimited memory at the moment of run
- Compiler: gcc (Spack GCC) 11.3.0
- Spack modules: cmake/3.26.3, gcc/11.3.0, intel-mkl/2019.3.199, intel-oneapi-mpi/2021.9.0
- GEOS: commit cb20b88c5e2d9afcf0b727bd8abb70f2c5366772
Screenshots
top (left) and bottom (right) views of properties
in MPa
in mD
If you are curious, here is the comparison of vertical displacements (subsidence) in poroelastic case at top layer. GEOS's displacemements are interpolated to cell-centers. 40x80x40 grid.
It seems there's a problem in the input file. I was expecting to see at least one WellControls block here: https://github.com/av-novikov/spe10_mech_geos/blob/2606b60196ba5473f88b7190a390005cd04d2c68/thermoporoelastic_fim.xml#L52-L56
@jhuang2601 can you help examining the input file?
Dear @victorapm,
Thanks for your hint about wells!
I added InternalWell in VTKMesh here: https://github.com/av-novikov/spe10_mech_geos/tree/with_internal_well.
However, still GMRES/MGR is stalling, UMFPACK runs out of memory.
To reproduce:
git clone -b with_internal_well https://github.com/av-novikov/spe10_mech_geos.gitcd spe10_mech_geos./path/to/geosx -i thermoporoelastic_fim.xml
Initialization of SinglePhasePoromechanicsReservoir works without wells on coarser mesh and in poroelastic setup.
However, if understand you correctly, non-empty targetRegions for well solver is a requirement?
Dear @jhuang2601,
Indeed <LinearSolverParameters directParallel="1"/> works either with or without well. I guess it takes reasonable, albeit quite some, time for 32k mesh. I have the following output
------------------- TIMESTEP START -------------------
- Time: -3168y, -319d, -4h-1m-4s (-100000000000 s)
- Delta Time: 3168y, 319d, 04h01m04s (100000000000 s)
- Cycle: 0
------------------------------------------------------
Time: -1.00e+11 s, dt: 100000000000 s, Cycle: 0
Task `singlePhasePoroelasticityPreEquilibrationStep`: at time -100000000000s, physics solver `singlePhasePoroelasticity` is set to perform stress initialization during the next time step(s)
Task `singlePhasePoroelasticityPreEquilibrationStep`: at time -100000000000s, physics solver `linearElasticity` is resetting total displacement and velocity to zero
Attempt: 0, ConfigurationIter: 0, NewtonIter: 0
( Rflow ) = ( 0.00e+00 ) ( Renergy ) = ( 0.00e+00 ) ( Rsolid ) = ( 1.72e+02 ) ( Rwell ) = ( 0.00e+00 ) ( R ) = ( 1.72e+02 )
Last LinSolve(iter,res) = ( 1, 4.48e-14 )
singlePhaseFlow: Max pressure change: 0.000 Pa (before scaling)
reservoirSystem: Global solution scaling factor = 1
Attempt: 0, ConfigurationIter: 0, NewtonIter: 1
( Rflow ) = ( 0.00e+00 ) ( Renergy ) = ( 0.00e+00 ) ( Rsolid ) = ( 1.73e-12 ) ( Rwell ) = ( 0.00e+00 ) ( R ) = ( 1.73e-12 )
wellControls1: well is shut
reservoirSystem: solver converged in less than 3 iterations, time-step required will be increased.
Not sure if possible to do the same for GMRES/MGR to calculate larger meshes? But thanks for a tip!
I see the issue now. We don't have an MGR strategy implemented for this physics model (thermalSinglePhasePoromechanicsReservoir). For now, one option is to impose source/sink fluxes via boxes and don't use wells. Another option is to keep the wells and run a compositional model instead of single-phase.
cc @castelletto1
Thanks @victorapm , compositional model works fine to me.
Not sure if absense of MGR strategy could explain the behavior of UMFPACK.
Let me know if you want to keep issue open due to MGR for thermalSinglePhasePoromechanicsReservoir.
Otherwise, i will close it.