Path: blob/devel/elmerice/examples/InverseMethods_OLD/Robin_Beta_GradientValid.sif
3203 views
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Compare the total derivative of the cost function computed as:
! (1) dJ=P.G where P is a perturbation vector of the variable of interest
! G is the gradient of the cost function computed by an inverse method
! (2) [J(V+hP)-J(V)]/h : forward finite difference computation of the derivative
! V is the variable of interest
! h is the step size
!
!
! Compute (1) from at the first iteration and update V=Vini+hP, h=1
! Compute (2) for all the other iteration with h^i+1=h^i/2
!
! RESULTS stored in :
! Gradient Validation section
! Result File = File "GradientValidation_$name".dat"
! gives: h, abs(1-2)/1, (1), (2)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
check keywords warn
! name of the run used for the outputs
$name="Robin_Beta_GradientValidation"
! Parametre de régularisation
$Lambda=0.0e00
! Domaine definition
$Lx = 200.0e3
$Ly = 50.0e03
!Reference Slip Coefficicient used to construct surface velocities
$ function betaSquare(tx) {\
Lx = 200.0e3;\
Ly = 50.0e03;\
yearinsec = 365.25*24*60*60;\
F1=sin(3.0*pi*tx(0)/Lx)*sin(pi*tx(1)/Ly);\
F2=sin(pi*tx(0)/(2.0*Lx))*cos(4.0*pi*tx(1)/Ly);\
beta=5.0e3*F1+5.0e03*F2;\
_betaSquare=beta*beta/(1.0e06*yearinsec);\
}
!Reference Viscosity used to construct surface velocities
$ function MuSquare(tx) {\
Lx = 200.0e3;\
Ly = 50.0e03;\
yearinsec = 365.25*24*60*60;\
F1=sin(3.0*pi*tx(0)/Lx)*sin(pi*tx(1)/Ly);\
F2=sin(pi*tx(0)/(2.0*Lx))*cos(4.0*pi*tx(1)/Ly);\
mu=sqrt(1.8e08)+1.0e03*(F1+F2);\
_MuSquare=mu*mu*1.0e-6*(2.0*yearinsec)^(-1.0/3.0);\
}
!some constants
$yearinsec = 365.25*24*60*60
$rhoi = 917.0/(1.0e6*yearinsec^2) ! MPa - a - m
$gravity = -9.81*yearinsec^2
Header
Mesh DB "." "mesh2D"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Simulation
Coordinate System = Cartesian 3D
Simulation Type = Steady State
! Internal extrusion
Extruded Mesh Levels=5
Output Intervals = 1
Steady State Max Iterations = 20
Steady State Min Iterations = 1
! Output File = "Test_$name".result"
! Post File = "Test_$name".ep"
Initialize Dirichlet Conditions = Logical False
max output level = 3
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Main ice body
Body 1
Equation = 1
Body Force = 1
Material = 1
Initial Condition = 1
End
! lower surface
Body 2
Equation = 2
Body Force = 1
Material = 1
Initial Condition = 1
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Initial Condition 1
BetaS = Variable coordinate 1, Coordinate 2
REAL MATC "betaSquare(tx)"
MuS = Variable coordinate 1, Coordinate 2
REAL MATC "MuSquare(tx)"
! initial guess for (square root) slip coeff.
Beta = REAL $ 1.0e3/sqrt(1.0e06*yearinsec)
! the perturbation
BetaP = REAL $ 1.0e3/sqrt(1.0e06*yearinsec)
Pressure = Real 0.0
Velocity 1 = Real 0.0
Velocity 2 = Real 0.0
Velocity 3 = Real 0.0
VeloD 1 = Real 0.0
VeloD 2 = Real 0.0
VeloD 3 = Real 0.0
VeloD 4 = Real 0.0
! Surface velocities (data)
Vsurfini 1 = Variable Coordinate 1
Real procedure "Executables/USF_Init" "UIni"
Vsurfini 2 = Variable Coordinate 1
Real procedure "Executables/USF_Init" "VIni"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body Force 1
Flow BodyForce 1 = Real 0.0
Flow BodyForce 2 = Real 0.0
Flow BodyForce 3 = Real $gravity
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! ice material properties in MPa - m - a system
Material 1
Density = Real $rhoi
Viscosity Model = String "power law"
Viscosity = Equals MuS
Viscosity Exponent = Real $1.0e00/3.0e00
Critical Shear Rate = Real 1.0e-10
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Solver 1
Exec Solver = Before simulation
Equation = "MapCoordinate"
Procedure = "StructuredMeshMapper" "StructuredMeshMapper"
Active Coordinate = Integer 3
End
!!!! Navier-Stokes Solution
Solver 2
Equation = "Navier-Stokes"
Stabilize = logical True
flow model = Stokes
!Linear System Solver = Direct
!Linear System Direct Method = mumps
!mumps percentage increase working space = integer 60
Linear System Solver = Iterative
Linear System Iterative Method = GMRES
Linear System GMRES Restart = 100
Linear System Preconditioning= ILU0
Linear System Convergence Tolerance= 1.0e-08
Linear System Max Iterations = 1000
!
Nonlinear System Max Iterations = Integer 100
Nonlinear System Convergence Tolerance = Real 1.0e-7
Nonlinear System Newton After Iterations = Integer 10
Nonlinear System Newton After Tolerance = Real 1.0e-03
Nonlinear System Relaxation Factor = Real 1.0
Nonlinear System Reset Newton = Logical True
Steady State Convergence Tolerance = Real 1.0e-12
! Define some useful Variables
Exported Variable 1 = BetaS
Exported Variable 1 DOFS = 1
Exported Variable 2 = MuS
Exported Variable 2 DOFS = 1
! square root of the slip coef
Exported Variable 3 = Beta
Exported Variable 3 DOFS = Integer 1
! perturbation variable
Exported Variable 4 = BetaP
Exported Variable 4 DOFS = Integer 1
! derivative of the cost fn wr to beta
Exported Variable 5 = DJDBeta
Exported Variable 5 DOFS = Integer 1
! value of the cost function
Exported Variable 6 = CostValue
Exported Variable 6 DOFS = Integer 1
Exported Variable 7 = VsurfIni
Exported Variable 7 DOFS = Integer 2
End
!!!! Navier-Stokes = Dirichlet Problem
Solver 3
Equation = "NS-Dirichlet"
Variable = VeloD
Variable Dofs = 4
procedure = "FlowSolve" "FlowSolver"
!Linear System Solver = Direct
!Linear System Direct Method = mumps
Linear System Solver = Iterative
Linear System Iterative Method = GMRES
Linear System GMRES Restart = 100
Linear System Preconditioning= ILU0
Linear System Convergence Tolerance= 1.0e-08
Linear System Max Iterations = 1000
Nonlinear System Max Iterations = Integer 100
Nonlinear System Convergence Tolerance = Real 1.0e-7
Nonlinear System Newton After Iterations = Integer 10
Nonlinear System Newton After Tolerance = Real 1.0e-03
Nonlinear System Relaxation Factor = Real 1.0
Nonlinear System Reset Newton = Logical True
Steady State Convergence Tolerance = Real 1.0e-12
End
!!! Compute Cost function
Solver 4
Equation = "Cost"
!! Solver need to be associated => Define dumy variable
Variable = -nooutput "CostV"
Variable DOFs = 1
procedure = "ElmerIceSolvers" "CostSolver_Robin"
Cost Variable Name = String "CostValue" ! Name of Cost Variable
Neumann Solution Name = String "Flow Solution"
Dirichlet Solution Name = String "VeloD"
Optimized Variable Name = String "Beta" ! Name of Beta for Regularization
Lambda = Real $Lambda ! Regularization Coef
! save the cost as a function of iterations
Cost Filename = File "Cost_$name".dat"
end
!!!!! Compute Derivative of Cost function / Beta
Solver 5
!!! needed only at the first iteration for computation of the total derivative with the inverse method
Exec Interval = 100
Equation = "DJDBeta"
!! Solver need to be associated => Define dumy variable
Variable = -nooutput "DJDB"
Variable DOFs = 1
procedure = "ElmerIceSolvers" "DJDBeta_Robin"
Neumann Solution Name = String "Flow Solution"
Dirichlet Solution Name = String "VeloD"
Optimized Variable Name = String "Beta" ! Name of Beta variable
Gradient Variable Name = String "DJDBeta" ! Name of gradient variable
PowerFormulation = Logical False
Beta2Formulation = Logical True ! SlipCoef define as Beta^2
Lambda = Real $Lambda ! Regularization Coef
end
!!!!! Gradient Validation
!!!!!! Compute total derivative and update the step size for the finite difference computation
Solver 6
Equation = "GradientValidation"
!! Solver need to be associated => Define dumy variable
Variable = -nooutput "UB"
Variable DOFs = 1
procedure = "./Executables/GradientValidation" "GradientValidation"
Cost Variable Name = String "CostValue"
Optimized Variable Name = String "Beta"
Perturbed Variable Name = String "BetaP"
Gradient Variable Name = String "DJDBeta"
Result File = File "GradientValidation_$name".dat"
end
Solver 7
Equation = "ResultOutput"
Procedure = File "ResultOutputSolve" "ResultOutputSolver"
Output File Name = string "Output_$name""
Vtu Format = logical true
Binary Output = True
Single Precision = True
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Equation 1
Active Solvers (5)= 1 2 3 4 7
NS Convect= False
End
Equation 2
Active Solvers (2)= 5 6
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Boundary Condition 1
Name = "Side Walls"
Target Boundaries(2) = 1 3
!Dirichlet BC
Velocity 1 = Real 0.0
Velocity 2 = Real 0.0
!Dirichlet BC => Same Dirichlet
VeloD 1 = Real 0.0
VeloD 2 = Real 0.0
End
Boundary Condition 2
Name = "Inflow"
Target Boundaries = 4
Velocity 1 = Variable Coordinate 2
REAL MATC "4.753e-6*yearinsec*(sin(2.0*pi*(Ly-tx)/Ly)+2.5*sin(pi*(Ly-tx)/Ly))"
Velocity 2 = Real 0.0
!Dirichlet BC => Same Dirichlet
VeloD 1 = Variable Coordinate 2
REAL MATC "4.753e-6*yearinsec*(sin(2.0*pi*(Ly-tx)/Ly)+2.5*sin(pi*(Ly-tx)/Ly))"
VeloD 2 = Real 0.0
End
Boundary Condition 3
Name = "Front"
Target Boundaries = 2
Velocity 1 = Variable Coordinate 2
REAL MATC "1.584e-5*yearinsec*(sin(2.0*pi*(Ly-tx)/Ly)+2.5*sin(pi*(Ly-tx)/Ly)+0.5*sin(3.0*pi*(Ly-tx)/Ly))"
Velocity 2 = Real 0.0
!Dirichlet BC => Same Dirichlet
VeloD 1 = Variable Coordinate 2
REAL MATC "1.584e-5*yearinsec*(sin(2.0*pi*(Ly-tx)/Ly)+2.5*sin(pi*(Ly-tx)/Ly)+0.5*sin(3.0*pi*(Ly-tx)/Ly))"
VeloD 2 = Real 0.0
End
Boundary Condition 4
!Name= "bed" mandatory to compute regularistaion term of the cost function (int (dbeta/dx) 2)
Name = "bed"
!Body Id used to solve
Body ID = Integer 2
Save Line = Logical True
Bottom Surface = Variable Coordinate 1
REAL procedure "Executables/USF_Init" "zbIni"
Normal-Tangential Velocity = Logical True
Normal-Tangential VeloD = Logical True
Velocity 1 = Real 0.0e0
VeloD 1 = Real 0.0e0
Slip Coefficient 2 = Variable Beta
REAL MATC "tx*tx"
Slip Coefficient 3 = Variable Beta
REAL MATC "tx*tx"
End
! Upper Surface
Boundary Condition 5
!Name= "Surface" mandatory to compute cost function
Name = "Surface"
Save Line = Logical True
! Used by StructuredMeshMapper for initial surface topography
! here interpolated from a regular DEM
Top Surface = Variable Coordinate 1
REAL procedure "Executables/USF_Init" "zsIni"
! Dirichlet problem applied observed velocities
VeloD 1 = Equals Vsurfini 1
VeloD 2 = Equals Vsurfini 2
End