Optimisation.sif

! Optimise the convection velocity to match the observed ie thicknesses
!-----------------------------
! PARAMETERS:
!----------------------------
$name=ipar(0)
$ID="Optim"

!# regularisation parameter
$LambdaReg=rpar(0)

!# H observations
$OBSERVATION_FILE="Hflight.txt"

!number of iterations and output intervals
$niter=200
$nout=25
!
include ../src/PARAMETERS.sif
!
Header
  Mesh DB "." "rectangle"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Simulation
  Coordinate System  = Cartesian 2D 

  Simulation Type = Steady

  Steady State Min Iterations = 1
  Steady State Max Iterations = $niter

  Post File = "$ID$_$name$.vtu"
  Output intervals = $nout

  Restart File = "Run0.result"
  Restart Before Initial Conditions = Logical True

  max output level = 3
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body 1
  Equation = 1
  Body Force = 1
  Initial Condition = 1
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! initial convection velocity
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Initial Condition 1
  Velocity 1 = Variable USolution 1, r 1
    REAL LUA "tx[0]+tx[1]"

  Velocity 2 = Variable USolution 2, r 2
    REAL LUA "tx[0]+tx[1]"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body Force 1
  Flow BodyForce 1 = Real 0.0                          
  Flow BodyForce 2 = Real 0.0
  Flow BodyForce 3 = Real $gravity

  Top Surface Accumulation = Equals SMB
  Bottom Surface Accumulation = Real 0.0
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! Compute H
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Solver 1
   Equation = "Thickness"
   Variable = -dofs 1 "H"

   Procedure = "ElmerIceSolvers" "AdjointThickness_ThicknessSolver"

   Linear System Solver = Direct
   Linear System Direct Method = umfpack

!!  the convection velocity (mean horizontal velocity)
    Flow Solution Name = String "Velocity"
!!    
   Exported Variable 1 = SMB
   Exported Variable 2 = HSolution
   Exported Variable 3 = -dofs 2 r
   Exported Variable 4 = -dofs 2 USolution
   Exported Variable 5 = -dofs 2 Velocity

   Exported Variable 6 = -global CostValue
   Exported Variable 7 = Hb

   Steady State Convergence Tolerance = Real 1.0e-12
End
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! Compute Cost
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Solver 2
  Equation = "Cost"
    procedure = "ElmerIceSolvers" "Adjoint_CostDiscSolver"

   Cost Variable Name = String "CostValue"  ! Name of Cost Variable
   Lambda = Real 1.0
 ! save the cost as a function of iterations (iterations,Cost,rms=sqrt(2*Cost/Ndata)
   Cost Filename = File "Cost_$name$.dat"

   Observed Variable Name = String "H"
   Observed Variable dimension = Integer 1

   Observation File Name = File "$OBSERVATION_FILE$"
end
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! Compute regularisation (penalise spatial derivatives of H)
!  has to be befoe the adjoint as in depends on H
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Solver 3
  Equation = String "CostReg"
  procedure = "ElmerIceSolvers" "Adjoint_CostRegSolver"

  !# True if cost function and gradient must be initialised to 0 in this solve
  Reset Cost Value = Logical False

  Cost Filename = File "CostReg_$name$.dat"

  Lambda=Real $LambdaReg

  Cost Variable Name = String "CostValue"
  Gradient Variable Name = String "Hb"
  Optimized Variable Name = String "H"

  A priori Regularisation = Logical False
end
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!  adjoint
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Solver 4
  Equation = "Adjoint"
  Variable = -nooutput Adjoint
  Variable Dofs = 1

  procedure = "ElmerIceSolvers" "Adjoint_LinearSolver"

!Name of the flow solution solver
  Direct Solver Equation Name = string "Thickness"

   Linear System Solver = Direct
   Linear System Direct Method = umfpack
End
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!  Gradient
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Solver 5
  Equation = "DJDp"

  procedure = "ElmerIceSolvers" "AdjointThickness_GradientSolver"

  Thickness Solution Name = String H
  Adjoint Solution Name = String "Adjoint"
  Flow Solution Name = String "Velocity"

  ComputeDJDUV = Logical True
  ComputeDJDsmbTop = Logical True

  Exported Variable 1 = DJDUV
  Exported Variable 1 DOFs = 2

  Exported Variable 2 = DJDsmbTop
  Exported Variable 2 DOFs = 1
end
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!  Optimisation
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Solver 6
  Equation = "Optimize_m1qn3"
  procedure = "ElmerIceSolvers" "Optimize_m1qn3Parallel"

  Cost Variable Name = String "CostValue"
  Optimized Variable Name = String "Velocity"
  Gradient Variable Name = String "DJDUV"
  gradient Norm File = String "GradientNormAdjoint_$name$.dat"

  Mesh Independent = Logical False

! M1QN3 Parameters
  M1QN3 dxmin = Real 1.0e-10
  M1QN3 epsg = Real  1.e-12
  M1QN3 niter = Integer $niter
  M1QN3 nsim = Integer $niter
  M1QN3 impres = Integer 5
  M1QN3 DIS Mode = Logical True
  M1QN3 df1 = Real 0.5
  M1QN3 normtype = String "dfn"
  M1QN3 OutputFile = File  "M1QN3_$name$.out"
  M1QN3 ndz = Integer 20
end
!!######################
Equation 1
  Active Solvers(6) = 1 2 3 4 5 6
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Boundary Condition 1
  name = "Sides"
  Target Boundaries(2) = 1 3

  ! we set dirichlet here as we may get inflow
  H = Equals HSolution
End

Boundary Condition 2
  name = "calving front"
  Target Boundaries = 2
  
End
Boundary Condition 3
  name = "inflow"
  Target Boundaries = 4

  H = Real $Hgl
End