Bedrock Sliding Boundary Condition

Extension of Elmer in computational glaciology
Post Reply
esvan
Posts: 3
Joined: 10 Aug 2015, 10:43
Antispam: Yes

Bedrock Sliding Boundary Condition

Post by esvan »

Hi, I’m a student and I’m approaching Elmer/Ice, trying to simulate the evolution of an existing glacier over time.
In particular, I use, as initial mesh, the extrusion in 10 layers of the glacier’s footprint, deformed following existing topographic data. As body force on the upper surface I give a summer ablation and a winter accumulation reading from a DAT file, and it works. A 0-velocity condition on the bedrock also works. But…

…I have a problem when I try to give a more realisic sliding condition on the bottom of the glacier. In particular with the help of Elmer courses’s material I tried to give a condition like this

Flow Force BC = Logical True
Normal-Tangential Velocity = Logical True
Velocity 1 = Real 0.0
Slip Coefficient 2 = Real 0.1
Slip Coefficient 3 = Real 0.1

But it just diverges. I just can’t solve this question. Anybody can kindly check my sif file to find anything I missed?

Here it is
sif_forum.txt
(8.1 KiB) Downloaded 603 times
Thank you
fgillet
Posts: 46
Joined: 30 Sep 2010, 16:58

Re: Bedrock Sliding Boundary Condition

Post by fgillet »

Hi,

Have you tried to use only Picard iterations (I see in your Navier-Stokes solver that you switch to Newton after 5 iterations)?
This is the first thing to try when the non linear iterations diverge.

cheers

Fabien
esvan
Posts: 3
Joined: 10 Aug 2015, 10:43
Antispam: Yes

Re: Bedrock Sliding Boundary Condition

Post by esvan »

Thank you for the reply. Following your suggestion I set the Navier-Stokes solver like following in order to avoid newton method

Nonlinear System Max Iterations = 50
Nonlinear System Convergence Tolerance = 1.0e-5
Nonlinear System Newton After Iterations = 50
Nonlinear System Newton After Tolerance = 1.0e-02
Nonlinear System Relaxation Factor = 1.00
Nonlinear System Reset Newton = Logical True

But unfortunately this didn't solve the problem. Looking at the terminal output I can see that the calculation diverges AFTER the first Navier-Stokes Solver. I post the terminal output below. It seems that it crashes when using my external function "getBalance". BUT that function works perfectly if I restore the boundary condition on the bedrock to "0-velocity". I attach also the code of my external functions named "utils.f90". I thank you and anybody else that can spend even a couple of minutes with my problem.

Here is the function code
utils.f90
(18.75 KiB) Downloaded 550 times
And here the output:

ELMER SOLVER (v 7.0) STARTED AT: 2015/08/22 11:53:45
MAIN:
MAIN: =============================================================
MAIN: ElmerSolver finite element software, Welcome!
MAIN: This program is free software licensed under (L)GPL
MAIN: Copyright 1st April 1995 - , CSC - IT Center for Science Ltd.
MAIN: Webpage http://www.csc.fi/elmer, Email elmeradm@csc.fi
MAIN: Library version: 7.0 (Rev: 1b38bc5)
MAIN: HYPRE library linked in.
MAIN: =============================================================
MAIN:
MAIN:
MAIN: -------------------------------------
MAIN: Reading Model: 1997-2007_test_Slip_22-08-2015.sif
Model Input: Unlisted keyword: [water density] in section: [constants]
Model Input: Unlisted keyword: [surfacedem] in section: [constants]
Model Input: Unlisted keyword: [bedrockdem] in section: [constants]
Model Input: Unlisted keyword: [switchperiod] in section: [constants]
Model Input: Unlisted keyword: [winterarrheniusfactor] in section: [constants]
Model Input: Unlisted keyword: [summerarrheniusfactor] in section: [constants]
Model Input: Unlisted keyword: [basalwintervelocity] in section: [constants]
Model Input: Unlisted keyword: [basalsummervelocity] in section: [constants]
Model Input: Unlisted keyword: [winteraccumulationfile] in section: [constants]
Model Input: Unlisted keyword: [summeraccumulationfile] in section: [constants]
Model Input: Unlisted keyword: [minicethickness] in section: [constants]
Model Input: Unlisted keyword: [zs accumulation] in section: [body force 2]
Loading user function library: [./utils_eros_fede]...[getBalance]
Model Input: Unlisted keyword: [zs] in section: [initial condition 2]
Loading user function library: [./utils_eros_fede]...[getSurfaceElevation]
Loading user function library: [./utils_eros_fede]...[getBedrockElevation]
Model Input: Unlisted keyword: [set arrhenius factor] in section: [material 1]
Model Input: Unlisted keyword: [arrhenius factor] in section: [material 1]
Loading user function library: [./utils_eros_fede]...[getArrheniusFactor]
Model Input: Unlisted keyword: [min zs] in section: [material 1]
Loading user function library: [./utils_eros_fede]...[getMinSurfaceElevation]
Model Input: Unlisted keyword: [max zs] in section: [material 1]
Model Input: Unlisted keyword: [mesh velocity variable] in section: [solver 1]
Model Input: Unlisted keyword: [mesh update variable] in section: [solver 1]
Model Input: Unlisted keyword: [mesh velocity first zero] in section: [solver 1]
Model Input: Unlisted keyword: [top surface variable name] in section: [solver 1]
Model Input: Unlisted keyword: [displacement mode] in section: [solver 1]
Model Input: Unlisted keyword: [correct surface] in section: [solver 1]
Model Input: Unlisted keyword: [minimum height] in section: [solver 1]
Model Input: Unlisted keyword: [stress variable name] in section: [solver 4]
Model Input: Unlisted keyword: [flow solver name] in section: [solver 4]
Model Input: Unlisted keyword: [relaxation factor] in section: [solver 6]
MAIN:
MAIN: -------------------------------------
MAIN: Time: 1/24 2592000.0000000000
MAIN: -------------------------------------
MAIN:
SolveEquations: -------------------------------------
SolveEquations: Coupled system iteration: 1
SolveEquations: -------------------------------------
StructuredMeshMapper: >Correct Surface< in case of intersecting upper and lower surface
StructuredMeshMapper: Adjusting upper surface to maintain minimum height to 0.5000E+00
getArrheniusFactor: Set value of Arrhenius Factor to 2.4000000E-24
ComputeChange: NS (ITER=1) (NRM,RELC): ( 0.66460863E-04 2.0000000 ) :: navier-stokes
ComputeChange: NS (ITER=2) (NRM,RELC): ( 0.10736614E-03 0.47064352 ) :: navier-stokes
ComputeChange: NS (ITER=3) (NRM,RELC): ( 0.14503338E-03 0.29847314 ) :: navier-stokes
ComputeChange: NS (ITER=4) (NRM,RELC): ( 0.17635048E-03 0.19488907 ) :: navier-stokes
ComputeChange: NS (ITER=5) (NRM,RELC): ( 0.20055361E-03 0.12843124 ) :: navier-stokes
ComputeChange: NS (ITER=6) (NRM,RELC): ( 0.21831422E-03 0.84802943E-01 ) :: navier-stokes
ComputeChange: NS (ITER=7) (NRM,RELC): ( 0.23089744E-03 0.56023585E-01 ) :: navier-stokes
ComputeChange: NS (ITER=8) (NRM,RELC): ( 0.23960695E-03 0.37021991E-01 ) :: navier-stokes
ComputeChange: NS (ITER=9) (NRM,RELC): ( 0.24554298E-03 0.24470931E-01 ) :: navier-stokes
ComputeChange: NS (ITER=10) (NRM,RELC): ( 0.24954789E-03 0.16178463E-01 ) :: navier-stokes
ComputeChange: NS (ITER=11) (NRM,RELC): ( 0.25223184E-03 0.10697718E-01 ) :: navier-stokes
ComputeChange: NS (ITER=12) (NRM,RELC): ( 0.25402263E-03 0.70747007E-02 ) :: navier-stokes
ComputeChange: NS (ITER=13) (NRM,RELC): ( 0.25726561E-03 0.12685526E-01 ) :: navier-stokes
ComputeChange: NS (ITER=14) (NRM,RELC): ( 0.25753787E-03 0.10577202E-02 ) :: navier-stokes
ComputeChange: NS (ITER=15) (NRM,RELC): ( 0.25755675E-03 0.73293926E-04 ) :: navier-stokes
ComputeChange: NS (ITER=16) (NRM,RELC): ( 0.25755835E-03 0.62131185E-05 ) :: navier-stokes
ComputeChange: SS (ITER=1) (NRM,RELC): ( 0.25755835E-03 2.0000000 ) :: navier-stokes
ComputeChange: NS (ITER=1) (NRM,RELC): ( 394880.79 2.0000000 ) :: sij
ComputeChange: NS (ITER=2) (NRM,RELC): ( 421559.88 0.65354617E-01 ) :: sij
ComputeChange: NS (ITER=3) (NRM,RELC): ( 431568.68 0.23463770E-01 ) :: sij
ComputeChange: NS (ITER=4) (NRM,RELC): ( 309568.48 0.32922437 ) :: sij
ComputeChange: NS (ITER=5) (NRM,RELC): ( 235878.41 0.27020070 ) :: sij
ComputeChange: NS (ITER=6) (NRM,RELC): ( 239879.73 0.16820801E-01 ) :: sij
ComputeChange: SS (ITER=1) (NRM,RELC): ( 239879.73 2.0000000 ) :: sij
ComputeChange: SS (ITER=1) (NRM,RELC): ( 0.0000000 0.0000000 ) :: eigenstresses
getBalance: Set accumulation file to balance_winter.dat
3 0.4927E-01
6 0.4202E+01
9 0.3234E+01
12 0.1999E+01
15 0.2004E+01
18 0.2006E+01
21 0.2006E+01
...
Cutted by me
...
1026 0.1892+139
1029 0.6784+140
1032 0.1688+142
1032 0.1688+142
ERROR:: IterSolve: Failed convergence tolerances.
fgillet
Posts: 46
Joined: 30 Sep 2010, 16:58

Re: Bedrock Sliding Boundary Condition

Post by fgillet »

Hi,

It's the iterative method for the linear system of the free surface (Solver 6) which is not converging.

Increasing the verbosity level (Max Output Level = I in the Simulation section; max I=42) can sometimes help to see where the problem comes from.
Have you checked your velocity field (you can set Exec Solver =never to your Free Surface solver to visualise your results of the NS solver)?
Units for the splip coeff in SI are Pa.s.m-1 (Tau=beta u). 0.1 in SI units is very low; you can strat by putting something like 1MPa.a.m-1 (be carefull my units), this will be very high friction (so you will get a solution close to your 0-velocity); then decrease the value.

In general:
- if iteratives methods for the linear system are not converging you can always switch to a direct solver to get a solution.
- Conditionning of the Free Surface solver is usually good and ietrative methods usually converge. I assume here that you have crazy velocities because of too low friction; and that why it is not converging
- The free surface equatiopn is linear if you are not using min/max (Apply dirichlet = False); so you need only 1 Non linear iteration to get the solution (Nonlinear System Max Iterations = 1, and Nonlinear System Relaxation Factor = 1).

Fabien
esvan
Posts: 3
Joined: 10 Aug 2015, 10:43
Antispam: Yes

Re: Bedrock Sliding Boundary Condition

Post by esvan »

Thank you so much for your answers.

We've solved the problem.
In particular I'm working in SI system, but the Slide Coefficient value of 0.1 was for MPa-m-a system. Calculating the appropriate value of coefficients gets a realistic velocity field and the Free Surface solver doesn't fail anymore.

Thank you again
Post Reply