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
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.