## Thermo-mechanically coupled Simulation question

Extension of Elmer in computational glaciology
jeffjoy
Posts: 23
Joined: 10 Nov 2020, 04:01
Antispam: Yes

### Thermo-mechanically coupled Simulation question

I have add the temperature solver,add the geothermal flux,however, the result .vtu have the different homolougous but have the same velocity,i think i don't Thermo-mechanically coupled successfully;but i don't know where have the mistakes;i use the exaples of 2021 elmerice courses synthetic glacier;the tempdiagnostic sif code as follows;who can help me,thanks?

Code: Select all

``````!---LUA BEGIN
! function IfThenElse(condition,t,f)
!   if condition then return t else return f end
! end
!---LUA END

!! in SI units, input in Kelvin
\$ function capacity(T) { _capacity=146.3+(7.253*T)}

!! in SI units, input in Kelvin
\$ function conductivity(T)  { _conductivity=9.828*exp(-5.7E-03*T)}

!! pressuremeltingpoint (Pressure in MPa)
\$ function pressuremeltingpoint(PIN) {\
P = PIN;\
if (P<0.0) P=0.0;\
beta=9.8E-08*1.0E06;\
_pressuremeltingpoint=273.15-(beta*P);\
}

\$namerun = "TempDiagnostic"
#slc=0.1
\$yearinsec = 365.25*24*60*60

Mesh DB "." "Mesh2d"
End

include "../Parameters/Physical_Parameters.IN"

!
!#Startyear = 0.0

!number of timesteps/year
#nt=10
!dt
#dt = 1.0/nt
!number of years
#ny=5

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Simulation
Coordinate System  =  Cartesian 3D

!Timestepping Method = "bdf"
!BDF Order = 2

!Timestep Intervals(1) = #nt*ny
!Output Intervals(1) = #5*nt
!Timestep Sizes(1) = #dt

Steady State Min Iterations = 1
Steady State Max Iterations = 50

!! Internal Extrusion
include Extrusion.sif

!! Restart
Restart Position = 0
!Restart Time = Real #Startyear
Restart Before Initial Conditions = Logical True

!! Output
Output File = "\$namerun\$.result"
Post File = "\$namerun\$.vtu"
vtu:VTU Time Collection = Logical True

!! scalar outputs to get timing infos
Scalars File = "Scalars_\$namerun\$.dat"
scalars: Parallel Reduce = Logical True
scalars: Variable 1 = String "Time"
scalars: Operator 2 = String "partitions"
scalars: Variable 3 = String "Flow Solution"
scalars: Operator 3 = String "dofs"

Simulation Timing = Logical True

max output level = 3
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body 1
Equation = 1
Body Force = 1
Material = 1
Initial Condition = 1
End
Body 2
Name= "surface"
Equation = 2
Material = 1
Body Force = 2
Initial Condition = 2
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Initial Condition 1
!Temperature = Real 271.0
End

Initial Condition 2
!Zs= Equals surfDEM
Ref Zs = Equals surfDEM

End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body Force 1
Flow BodyForce 1 = Real 0.0
Flow BodyForce 2 = Real 0.0
Flow BodyForce 3 = Real #gravity

Temperature Volume Source = Equals W ! The volumetric heat source

!Temperature Passive = Variable "Thickness","slope"
!  Real LUA "IfThenElse((tx[0] <= 50) or (tx[1] > 1.2), 1.0, -1.0)"

! This should be in Body Force 2 but not working
! for solver executed on a boundary
Zs = Variable bedDEM
Real LUA "tx[0]+ MinH"
! should make it also dependant of SMB, i.E. allow zs to change if SMB>0
Zs Condition = Variable "IcyMask","Mass Balance"
Real LUA "IfThenElse((tx[0]< -0.5) and (tx[1] <= 0.0), 1.0, -1.0)"
End

Body Force 2
Zs Accumulation Flux 1 = real 0.0
Zs Accumulation Flux 2 = real 0.0
Zs Accumulation Flux 3 = Equals "Mass Balance"

! surface slope norm
slope = Variable "dzs 1", "dzs 2"
REAL LUA "math.sqrt(tx[0]*tx[0]+tx[1]*tx[1])"

! mask mass balance with surface slope
Mass Balance = Variable "Mass Balance Ini","slope"
Real LUA "IfThenElse(tx[0]>0,IfThenElse(tx[1]< 1.2, tx[0], 0.0),tx[0])"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Material 1
! For the ice flow
Density = Real #rhoi

Viscosity Model = String "Glen"
Viscosity = Real 1.0
Glen Exponent = Real 3.0
Critical Shear Rate = Real 1.0e-10

! properties with T
Rate Factor 1 = Real #A1
Rate Factor 2 = Real #A2
Activation Energy 1 = Real #Q1
Activation Energy 2 = Real #Q2
Glen Enhancement Factor = Real 1.0
Limit Temperature = Real -10.0
!Relative Temperature = Real 0.0

! or provide A
!Set Arrhenius Factor = Logical True
!Arrhenius Factor = Real #A

! the heat capacity as a MATC function of temperature itself
!-----------------------------------------------------------
Temperature Heat Capacity = Variable Temperature
Real MATC "capacity(tx)*yearinsec^2"
! the heat conductivity as a MATC function of temperature itself
!--------------------------------------------------------------
Temperature Heat Conductivity = Variable Temperature
Real MATC "conductivity(tx)*yearinsec*1.0E-06"
! Upper limit - pressure melting point
!  as a MATC function of the pressure (what else?)
!-------------------------------------------------
Temperature Upper Limit = Variable Pressure
Real MATC "pressuremeltingpoint(tx)"
! lower limit (to be save) as 0 K
!--------------------------------
Temperature Lower Limit = Real 0.0

Min Zs = Variable bedDEM
Real LUA "tx[0]+ MinH"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Map mesh to bottom and top surfaces
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Solver 1
Equation = "MapCoordinate"
Procedure = "StructuredMeshMapper" "StructuredMeshMapper"

! Extrusion direction
Active Coordinate = Integer 3

Displacement Mode = Logical False

! Check for critical thickness
Correct Surface = Logical True
Minimum Height = Real #MinH

! Top and bottom surfaces defined from variables
Top Surface Variable Name = String "Zs"
Bottom Surface Variable Name = String "bedDEM"
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute Thickness and Depth
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Solver 2
Equation = "HeightDepth 1"
Procedure = "StructuredProjectToPlane" "StructuredProjectToPlane"

Active Coordinate = Integer 3

Project to everywhere = Logical True

Operator 1 = Thickness
Operator 2 = Depth
!  Operator 3 = Height
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Solver 3
! to be executed on top surface (need Thickness)
Variable DOFs = 1

! no matrix resolution
Optimize Bandwidth = Logical False

Toler = Real 1.0e-1
Ice Free Thickness = Real #MinH
Remove Isolated Points = Logical True
Remove Isolated Edges = Logical True
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Stokes
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Solver 4
Equation = "Stokes-Vec"
Procedure = "IncompressibleNSVec" "IncompressibleNSSolver"
Stokes Flow = logical true
!Exec Solver = "Never"

Div-Curl Discretization = Logical False
Stabilization Method = String Stabilized

!linear settings:
!------------------------------
include linsys/BiCGStab.sif

!Non-linear iteration settings:
!------------------------------
Nonlinear System Max Iterations = 20
Nonlinear System Convergence Tolerance  = 1.0e-5
Nonlinear System Newton After Iterations = 4
Nonlinear System Newton After Tolerance = 1.0e-4
Nonlinear System Reset Newton = Logical True
Nonlinear System Abort Not Converged = Logical True

! Convergence on timelevel (not required here)
!---------------------------------------------
Steady State Convergence Tolerance = Real 1.0e-3

Relative Integration Order = -1
!Number of Integration Points = Integer 44 ! 21, 28, 44, 64, ...

! 1st iteration viscosity is constant
Constant-Viscosity Start = Logical False

!! Timing infos
Bulk Assembly Timing = Logical True
Linear System Timing = Logical True
Linear System Timing Cumulative = Logical True
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Temperature
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Solver 5
Equation = DeformationalHeat
Variable = W
Variable DOFs = 1

procedure =  "ElmerIceSolvers" "DeformationalHeatSolver"

Linear System Solver = direct
Linear System direct Method = umfpack
End

Solver 6
Equation = String "Homologous Temperature Equation"
Procedure =  File "ElmerIceSolvers" "TemperateIceSolver"
! Comment next line in parallel, as EliminateDirichlet does
! not work in parallel
!------------------------------------------------------------
Loop While Unconstrained Nodes = Logical True
Variable = String "Temperature"
Variable DOFs = 1
Linear System Solver = "Iterative"
Linear System Iterative Method = "BiCGStab"
Linear System Max Iterations = 1000
Linear System Convergence Tolerance = 1.0E-07
Linear System Abort Not Converged = True
Linear System Preconditioning = "ILU6"
Linear System Residual Output = 1
Steady State Convergence Tolerance = 1.0E-04
Nonlinear System Convergence Tolerance = 1.0E-05
Nonlinear System Max Iterations = 50
Nonlinear System Relaxation Factor = Real 9.999E-01
! uses the contact algorithm (aka Dirichlet algorithm)
!-----------------------------------------------------
Apply Dirichlet = Logical True
Stabilize = True
! those two variables are needed in order to store
! the relative or homologous temperature as well
! as the residual
!-------------------------------------------------
Exported Variable 1 = String "Temperature Homologous"
Exported Variable 1 DOFs = 1
Exported Variable 2 = String "Temperature Residual"
Exported Variable 2 DOFs = 1
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Surface slope
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Solver 7
Equation = "Surface Slope"
Variable = -nooutput "dzs"
Variable DOFs = 2

Variable Name = String "zs"

Update Exported Variables = Logical True
Exported Variable 1 = "slope"
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Mass balance forcing
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Solver 8
Equation = SMB
Procedure = "../SRC/SyntSMB" "SyntSMB"
Variable = "Mass Balance Ini"
Variable DOFs = 1
Exec Solver = Never

! no matrix resolution
Optimize Bandwidth = Logical False

SMB Exponent = Real #SMBExponent
SMB ELA = Real #SMBELA

Update Exported Variables = Logical True
Exported Variable 1 = "Mass Balance"
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Free surface evolution
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Solver 9
Equation =  String "Free Surface Evolution"
Procedure = "FreeSurfaceSolver" "FreeSurfaceSolver"
Exec Solver = Never

Variable = "Zs"
Variable DOFs = 1

! calculate dz/dt (better than from mesh velocity in case os steady-state iterations)
Calculate Velocity = Logical True

! Apply internal limiters
Apply Dirichlet = Logical true

! Steb method
Stabilization Method = Stabilized

! linear settings
Linear System Solver = Iterative
Linear System Iterative Method = BiCGStab
Linear System Max Iterations  = 1000
Linear System Preconditioning = ILU0
Linear System Convergence Tolerance = 1.0e-10

! non-linear settings
Nonlinear System Max Iterations = 20 ! variational inequality needs more than one round
Nonlinear System Min Iterations = 2
Nonlinear System Convergence Tolerance = 1.0e-8

Steady State Convergence Tolerance = 1.0e-6

! loads also takes into account dirichlet conditions
! to compute residual flux

Exported Variable 1 = -nooutput "Zs Residual"
Exported Variable 2 =  "Ref Zs"
End

Solver 10
! to be executed on top surface (need Thickness)
Exec Solver = After Timestep
Equation = "Save 1D Vars"
Procedure = "ElmerIceSolvers" "Scalar_OUTPUT"
Variable = -nooutput "savescal"

File Name = File "1DVar_OUTPUT_\$namerun\$.dat"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Equation 1
Active Solvers(5) = 1 2 4 5 6
Flow Solution Name = String "Flow Solution"
Convection = Computed
End

Equation 2
Active Solvers(5) = 3 7 8 9 10
Flow Solution Name = String "Flow Solution"
Convection = Computed
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Side
Boundary Condition 1
Target Boundaries = 1
Name = "side"

Normal-Tangential Velocity = Logical True
Velocity 1 = Real 0.0e0

End

! Bedrock
Boundary Condition 2
Name = "bed"

!-------------------
! geothermal heatflux
!--------------------
Temperature Flux BC = Logical True
Temperature Heat Flux = Real \$56.05E-03*365.25*24.0*3600.0*1.0E-6
!-------------------
! frictional heat
!--------------------
Temperature Load = Variable Velocity 1

include BCs/noslip.sif

End

! Upper Surface
Boundary Condition 3
Name = "upper surface"
Body Id = 2

Temperature = Real 273.0
End

``````
jeffjoy
Posts: 23
Joined: 10 Nov 2020, 04:01
Antispam: Yes

### Re: Thermo-mechanically coupled Simulation question

Nobody knows? Is there a bug of the ELmer/ice?
jeffjoy
Posts: 23
Joined: 10 Nov 2020, 04:01
Antispam: Yes

### Re: Thermo-mechanically coupled Simulation question

If i can't solve the problem ,my team will change to another model
Rich_B
Posts: 407
Joined: 24 Aug 2009, 20:18

### Re: Thermo-mechanically coupled Simulation question

Please post an archive of a working (or not working) solution.

Include the sif file and enough of the geometry input files to allow someone to fully run your simulation.
If possible, it also helps to include a copy of the elmersolver output in the archive.