Path: blob/devel/elmerice/examples/InverseMethods_OLD/Adjoint_Mu_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="Adjoint_Mu_GradientValidation"
$Lambda=0.0e00 ! Regularistaion parameter! not implemented for Mu optimisation!! has to be 0!!
! 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
Constants
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"
max output level = 3
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Main ice body
Body 1
Equation = 1
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) the ice viscosity
Mu = REAL $sqrt(1.8e8*1.0e-6*(2.0*yearinsec)^(-1.0/3.0))
! the perturbation variable
MuP = REAL $0.1*sqrt(1.8e8*1.0e-6*(2.0*yearinsec)^(-1.0/3.0))
Pressure = Real 0.0
Velocity 1 = Real 0.0
Velocity 2 = Real 0.0
Velocity 3 = Real 0.0
Adjoint 1 = Real 0.0
Adjoint 2 = Real 0.0
Adjoint 3 = Real 0.0
Adjoint 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 9.1376e-19
Viscosity Model = String "power law"
! Viscosity defined as mu^2 to ensure >0
Viscosity = Variable mu
Real MATC "tx*tx"
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 = forward model
Solver 2
Equation = "Navier-Stokes"
Stabilize = logical True
flow model = Stokes
!mandatory to save bulk stiffness matrix
calculate loads = Logical True
Linear System Solver = Iterative
Linear System Iterative Method = GMRES
Linear System GMRES Restart = 100
Linear System Preconditioning= ILU0
Linear System Convergence Tolerance= 1.0e-12
Linear System Max Iterations = 1000
! system self adjoint if Newton is used for the last iterations
Nonlinear System Max Iterations = Integer 100
Nonlinear System Convergence Tolerance = Real 1.0e-10
Nonlinear System Newton After Iterations = Integer 10
Nonlinear System Newton After Tolerance = Real 1.0e-03
Nonlinear System Relaxation Factor = Real 1.0
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 Viscosity
Exported Variable 3 = mu
Exported Variable 3 DOFS = Integer 1
! the perturbation variable
Exported Variable 4 = muP
Exported Variable 4 DOFS = Integer 1
! derivative of the cost fn wr to mu
Exported Variable 5 = DJDmu
Exported Variable 5 DOFS = Integer 1
! value of the cost function
Exported Variable 6 = CostValue
Exported Variable 6 DOFS = Integer 1
! Variable used to force the adjoint system/update in CostSolver
Exported Variable 7 = Velocityb
Exported Variable 7 DOFs = Integer 4
Exported Variable 8 = VsurfIni
Exported Variable 8 DOFS = Integer 2
End
!!! Compute Cost function
!!!!!!!! Has to be run before the Adjoint Solver as adjoint forcing is computed here !!!!!
Solver 3
Equation = "Cost"
!! Solver need to be associated => Define dumy variable
Variable = -nooutput "CostV"
Variable DOFs = 1
procedure = "ElmerIceSolvers" "CostSolver_Adjoint"
Cost Variable Name = String "CostValue" ! Name of Cost Variable
Optimized Variable Name = String "Mu" ! 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
!!!! Adjoint Solver
Solver 4
!!! needed only at the first iteration for computation of the total derivative with the inverse method
Exec Interval = 100
Equation = "Adjoint"
Variable = Adjoint
Variable Dofs = 4
procedure = "ElmerIceSolvers" "AdjointSolver"
!Name of the flow solution solver
Flow Solution Equation Name = string "Navier-Stokes"
Linear System Solver = Iterative
Linear System Iterative Method = GMRES
Linear System GMRES Restart = 100
Linear System Preconditioning= ILU0
Linear System Convergence Tolerance= 1.0e-12
Linear System Max Iterations = 1000
End
!!!!! Compute Derivative of Cost function / Mu
Solver 5
!!! needed only at the first iteration for computation of the total derivative with the inverse method
Exec Interval = 100
Equation = "DJDmu"
!! Solver need to be associated => Define dumy variable
Variable = -nooutput "DJDBV"
Variable DOFs = 1
procedure = "DJDmu_Adjoint" "DJDMu_Adjoint"
Flow Solution Name = String "Flow Solution"
Adjoint Solution Name = String "Adjoint"
Optimized Variable Name = String "mu" ! Name of Beta variable
Gradient Variable Name = String "DJDmu" ! Name of gradient variable
SquareFormulation = Logical True ! Viscosity not define as Mu^2
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 "mu"
Perturbed Variable Name = String "muP"
Gradient Variable Name = String "DJDmu"
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
Vector Field 1 = "Velocity"
Scalar Field 1 = Mu
Scalar Field 2 = MuS
Scalar Field 3 = Adjoint 1
Scalar Field 4 = Adjoint 2
Scalar Field 5 = Adjoint 3
Scalar Field 6 = VsurfIni 1
Scalar Field 7 = VsurfIni 2
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Equation 1
Active Solvers (7)= 1 2 3 4 5 6 7
NS Convect= False
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Boundary Condition 1
Name = "Side Walls"
Target Boundaries(2) = 1 3
Velocity 1 = Real 0.0
Velocity 2 = Real 0.0
!Dirichlet BC => Dirichlet = 0 for Adjoint
Adjoint 1 = Real 0.0
Adjoint 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 => Dirichlet = 0 for Adjoint
Adjoint 1 = Real 0.0
Adjoint 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 => Dirichlet = 0 for Adjoint
Adjoint 1 = Real 0.0
Adjoint 2 = Real 0.0
End
Boundary Condition 4
Name = "bed"
!Target Boundaries = 5
Save Line = Logical True
Bottom Surface = Variable Coordinate 1
REAL procedure "Executables/USF_Init" "zbIni"
Normal-Tangential Velocity = Logical True
Normal-Tangential Adjoint = Logical True
Adjoint Force BC = Logical True
Velocity 1 = Real 0.0e0
Adjoint 1 = Real 0.0e0
Slip Coefficient 2 = Equals BetaS
Slip Coefficient 3 = Equals BetaS
End
! Upper Surface
Boundary Condition 5
Name = "Surface"
!Target Boundaries = 6
Save Line = Logical True
Top Surface = Variable Coordinate 1
REAL procedure "Executables/USF_Init" "zsIni"
! Definition of the Cost function
Adjoint Cost = Variable Velocity 1 , Vsurfini 1 , Velocity 2 , Vsurfini 2
Real MATC "0.5*((tx(0)-tx(1))*(tx(0)-tx(1))+(tx(2)-tx(3))*(tx(2)-tx(3)))"
! derivative of the cost function wr u and v
Adjoint Cost der 1 = Variable Velocity 1 , Vsurfini 1
Real MATC "tx(0)-tx(1)"
Adjoint Cost der 2 = Variable Velocity 2 , Vsurfini 2
Real MATC "tx(0)-tx(1)"
End