Merge branch 'elmerice' into devel

elmerice/IceSheet/Greenland/SSA_FRICTION_INIT/OPTIM_BETA.sif

@@ -61,19 +61,15 @@ Body Force 1
   Flow BodyForce 2 = Real 0.0              
   Flow BodyForce 3 = Real $gravity
 
-! change of variable; DJDBeta is the derivative of the Cost fn w.r.t. the slip coeff.
-  DJDalpha = Variable DJDBeta , alpha
-     REAL procedure "ElmerIceUSF" "Derivative_TenPowerA"
-
 !# For Jdiv
   Top Surface Accumulation = Equals smb
   Bottom Surface Accumulation = Real 0.0
   Observed dhdt = Equals dhdt_obs
 
 !# Cost not computed if H<=Hmin
-  CostV Passive = Variable H
+  Cost_var Passive = Variable H
          Real procedure "USFs" "PassiveCond_H"
-  CostDiv Passive = Variable H
+  Cost_Cost_DHDT Passive = Variable H
          Real procedure "USFs" "PassiveCond_H"
   Passive Element Min Nodes = Integer 3
 
@@ -95,6 +91,8 @@ Material 1
   ! The friction parameter is 10^the optimised variable to insure > 0
   SSA Friction Parameter = Variable alpha
       REAL procedure "ElmerIceUSF" "TenPowerA"
+  SSA Friction Parameter Derivative = Variable alpha
+     REAL procedure "ElmerIceUSF" "TenPowerA_d"
 End
 !#######################################################
 !#######################################################
@@ -104,18 +102,17 @@ Solver 1
 
   Procedure = "ElmerIceSolvers" "AdjointSSA_SSASolver"
 
-!! Mandatory for the adjoint 
-  Calculate Loads = Logical True
-
 !! NUMERICAL INFORMATION
   !Linear System Solver = Direct
   !Linear System Direct Method = mumps
+  !mumps percentage increase working space = integer 40
    Linear System Solver = Iterative
-   Linear System Max Iterations = 300
-   Linear System Iterative Method = BiCGStab
+   Linear System Max Iterations = 500
+   Linear System Iterative Method = GCR
+   Linear System GCR Restart =  250
    Linear System Preconditioning = ILU2
    Linear System Abort Not Converged = False
-   Linear System Residual Output = 1000
+   Linear System Residual Output = 500
    Linear System Convergence Tolerance = 1.0e-12
 
    Nonlinear System Max Iterations = 30
@@ -131,56 +128,39 @@ Solver 1
   GL integration points number = Integer 20
 !!!!!
 !! Variables required for the inverse method
-  Exported Variable 1 = -nooutput CostValue
-  Exported Variable 2 = -nooutput DJDBeta
-  Exported Variable 3 = -nooutput "Velocityb"
-  Exported Variable 3 DOFs = 2
-  Exported Variable 4 = -nooutput "ssavelocity loads"
+  Exported Variable 1 = -global CostValue
+  Exported Variable 2 = -nooutput alpha
+  Exported Variable 3 = -nooutput DJDalpha
+  Exported Variable 4 = -nooutput "Velocityb"
   Exported Variable 4 DOFs = 2
 End
 !#######################################################
 !!! Compute Cost function
 !! Here the cost is the discrete sum_1^Ndata 1/2 ||u-u^obs|| evaluated at the data location (which may not correspond to mesh nodes)
 Solver 2
   Equation = "Cost"
-  !!  Solver need to be associated => Define dumy variable
-    Variable = -nooutput "CostV"
-    Variable DOFs = 1
-
-    procedure = "ElmerIceSolvers" "AdjointSSA_CostDiscSolver"
+  procedure = "ElmerIceSolvers" "Adjoint_CostDiscSolver"
 
-   Optimize Bandwidth = logical false
-
-   Problem Dimension = Integer 2 !2D mesh and 2D SSA Solution
    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_U_$name$.dat"
-   Lambda=Real 1.0
+   Cost Filename = File "Cost_$name$.dat"
 
    Observed Variable Name = String "SSAVelocity"
- ! File with data: 
- !   ASCII x,y,u,v
- ! or netcdf with vx, vy
-   Observation File Name = File "$VELOCITY_DATA$"
- ! 
-   Save used data = logical False
+   Observed Variable dimension = Integer 2
+
+ ! ASCII File with data: x,y,u,v
+   Observation File Name = File "$OBSERVATION_FILE$"
 end
 !#######################################################
 !!! Compute Cost function Jdiv (mismatch between model and observed flux divergence)
 !!  Jdiv=int_ 0.5 * (dhdt_obs + div(uH)-MB)^2
 Solver 3
   Equation = "Cost_DHDT"
-  !!  Solver need to be associated => Define dumy variable
-    Variable = -nooutput "CostDiv"
-    Variable DOFs = 1
-
    procedure = "ElmerIceSolvers" "AdjointSSA_CostFluxDivSolver"
 
-   Optimize Bandwidth = logical false
-
    Reset Cost Value = Logical False
 
-   Problem Dimension = Integer 2 !2D mesh and 2D SSA Solution
    Cost Variable Name = String "CostValue"  ! Name of Cost Variable
 
    Lambda= Real $LambdaDiv
@@ -198,106 +178,85 @@ Solver 4
   Variable = -nooutput Adjoint
   Variable Dofs = 2
 
-  procedure = "ElmerIceSolvers" "AdjointSSA_AdjointSolver"
+  procedure = "ElmerIceSolvers" "Adjoint_LinearSolver"
 
 !Name of the flow solution solver
-  Flow Solution Equation Name = string "SSA"
+   Direct Solver Equation Name = string "SSA"
+
+!  Linear System Solver = Direct
+!  Linear System Direct Method = mumps
+!  mumps percentage increase working space = integer 40
+
+  Linear System Solver = Iterative
+  Linear System Max Iterations = 500
+  Linear System Iterative Method = GCR
+  Linear System GCR Restart =  250
+  Linear System Preconditioning = ILU2
+  Linear System Abort Not Converged = False
+  Linear System Residual Output = 500
+  Linear System Convergence Tolerance = 1.0e-12
 
-!! Numerical Stuff
-  !Linear System Solver = Direct
-  !Linear System Direct Method = mumps
-   Linear System Solver = Iterative
-   Linear System Max Iterations = 300
-   Linear System Iterative Method = BiCGStab
-   Linear System Preconditioning = ILU2
-   Linear System Abort Not Converged = False
-   Linear System Residual Output = 1000
-   Linear System Convergence Tolerance = 1.0e-12
 End
 !#######################################################
 !!!!!  Compute Derivative of Cost function / Beta
 Solver 5
   Equation = "DJDBeta"
-
-  !!  Solver need to be associated => Define dumy variable
-    Variable = -nooutput "DJDB"
-    Variable DOFs = 1
-
-    procedure = "ElmerIceSolvers" "AdjointSSA_GradientSolver"
-
-    Optimize Bandwidth = logical false
+   procedure = "ElmerIceSolvers" "AdjointSSA_GradientSolver"
 
     Flow Solution Name = String "SSAVelocity"
     Adjoint Solution Name = String "Adjoint"
-
-    Compute DJDBeta = Logical True   ! Derivative with respect to the Friction parameter
+    ! Derivative with respect to the Friction parameter
+    ! here will be with respect to alpha (see Material)
+    Compute DJDBeta = Logical True
+    DJDBeta Name = String "DJDalpha"
     Reset DJDBeta = Logical True
-
 end
 !#######################################################
-Solver 6
-  Equation = "UpdateExport"
-   Variable = -nooutput "dumy"
-    Procedure = File "ElmerIceSolvers" "UpdateExport"
-    Optimize Bandwidth = logical false
-
-!used here to update DJDalpha from DJDbeta (see corresponding line in Body Force section)
-  Exported Variable 1 =  alpha
-  Exported Variable 2 =  -nooutput DJDalpha
-End
-!#######################################################
-!!!!!  Compute Regularistaion term
+!!!!!  Compute Regularisation term
 !   Regularisation by default is: Lambda * int_{Pb dimension} 0.5 * (d(var)/dx)**2 
 !     OUTPUT are : J and DJDvar
-Solver 7
-  Equation = "DJDBeta_Reg"
-
-  !!  Solver need to be associated => Define dumy variable
-    Variable = -nooutput "DJDBReg"
-    Variable DOFs = 1
-
-    procedure = "ElmerIceSolvers" "AdjointSSA_CostRegSolver"
-    Optimize Bandwidth = logical false
+Solver 6
+  Equation = "CostReg"
+   procedure = "ElmerIceSolvers" "Adjoint_CostRegSolver"
 
-    Problem Dimension=Integer 2
     Cost Filename=File "CostReg_$name$.dat"
     Optimized Variable Name= String "alpha"
     Gradient Variable Name= String "DJDalpha"
     Cost Variable Name= String "CostValue"
     Lambda= Real $LambdaReg
     Reset Cost Value= Logical False  !=> DJDapha already initialized in solver DJDBeta; switch off initialisation to 0 at the beginning of this solver
     A priori Regularisation= Logical False
+
 end
 !#######################################################
-!!!!! Optimization procedure : Parallel only
-Solver 8
+!!!!! Optimization procedure :
+Solver 7
   Equation = "Optimize_m1qn3"
-  !!  Solver need to be associated => Define dumy variable
-    Variable =  -nooutput "UB"
-    Variable DOFs = 1
-
   procedure = "ElmerIceSolvers" "Optimize_m1qn3Parallel"
 
   Cost Variable Name = String "CostValue"
   Optimized Variable Name = String "alpha"
   Gradient Variable Name = String "DJDalpha"
-  gradient Norm File = String "GradientNormAdjoint_$name$.dat"
+  gradient Norm File = File "GradientNormAdjoint_$name$.dat"
+
+  !!
+  Mesh Independent = Logical FALSE
 
  ! M1QN3 Parameters
   M1QN3 dxmin = Real 1.0e-10
-  M1QN3 epsg = Real  1.e-8
+  M1QN3 epsg = Real  1.e-6
   M1QN3 niter = Integer $niter
   M1QN3 nsim = Integer $niter
   M1QN3 impres = Integer 5
-  M1QN3 DIS Mode = Logical False
-  M1QN3 df1 = Real 0.1
+  M1QN3 DIS Mode = Logical true
+  M1QN3 df1 = Real 0.5
   M1QN3 normtype = String "dfn"
-  M1QN3 OutputFile = File  "M1QN3_$name$.dat"
+  M1QN3 OutputFile = File  "M1QN3_$name$.out"
   M1QN3 ndz = Integer 20
 
 end
 !#######################################################
-Solver 9
+Solver 8
   Equation = "UpdateExport2"
    Variable = -nooutput "dumy2"
     Procedure = File "ElmerIceSolvers" "UpdateExport"
@@ -308,7 +267,7 @@ End
 !#######################################################
 !#######################################################
 Equation 1
-  Active Solvers(9) = 1 2 3 4 5 6 7 8 9
+  Active Solvers(8) = 1 2 3 4 5 6 7 8
 End
 !#######################################################
 !#######################################################

elmerice/Solvers/AdjointSSA/AdjointSSA_AdjointSolver.F90

@@ -104,6 +104,15 @@ SUBROUTINE AdjointSSA_AdjointSolver( Model,Solver,dt,TransientSimulation )
    LOGICAL, SAVE :: Firsttime=.TRUE.
 
 
+   CALL Info(SolverName,'***********************',level=0)
+   CALL Info(SolverName,' This solver has been replaced by:',level=0)
+   CALL Info(SolverName,'   Adjoint_LinearSolver  ',level=0)
+   CALL Info(SolverName,' See documentation under:   ',level=0)
+   CALL Info(SolverName,'   elmerice/Solvers/Documentation   ',level=0)
+   CALL Info(SolverName,'***********************',level=0)
+   CALL FATAL(SolverName,' Use new solver !!')
+
+
    DIM = CoordinateSystemDimension()
 
    StiffMatrix => Solver % Matrix

elmerice/Solvers/AdjointSSA/AdjointSSA_CostContSolver.F90

@@ -70,7 +70,8 @@ SUBROUTINE AdjointSSA_CostContSolver( Model,Solver,dt,TransientSimulation )
   LOGICAL :: TransientSimulation
 !  
   CHARACTER(LEN=MAX_NAME_LEN), PARAMETER :: DefaultCostFile = 'CostOfT.dat'
-  CHARACTER(LEN=MAX_NAME_LEN) :: SolverName,CostFile
+  CHARACTER(LEN=MAX_NAME_LEN) :: SolverName='CostSolver_Adjoint'
+  CHARACTER(LEN=MAX_NAME_LEN) :: CostFile
   CHARACTER(LEN=MAX_NAME_LEN) :: CostSolName
   TYPE(Element_t),POINTER ::  Element
   TYPE(Variable_t), POINTER :: TimeVar,CostVar
@@ -94,6 +95,13 @@ SUBROUTINE AdjointSSA_CostContSolver( Model,Solver,dt,TransientSimulation )
   save SolverName,CostSolName,CostFile
   save ElementNodes
 
+   CALL Info(SolverName,'***********************',level=0)
+   CALL Info(SolverName,' This solver has been replaced by:',level=0)
+   CALL Info(SolverName,'   Adjoint_CostContSolver  ',level=0)
+   CALL Info(SolverName,' See documentation under:   ',level=0)
+   CALL Info(SolverName,'   elmerice/Solvers/Documentation   ',level=0)
+   CALL Info(SolverName,'***********************',level=0)
+   CALL FATAL(SolverName,' Use new solver !!')
 
   SolverParams => GetSolverParams()
   DIM=GetInteger(SolverParams ,'Problem Dimension',Found)

elmerice/Solvers/AdjointSSA/AdjointSSA_CostDiscSolver.F90

@@ -75,7 +75,8 @@ SUBROUTINE AdjointSSA_CostDiscSolver( Model,Solver,dt,TransientSimulation )
   LOGICAL :: TransientSimulation
 !  
   CHARACTER(LEN=MAX_NAME_LEN), PARAMETER :: DefaultCostFile = 'CostOfT.dat'
-  CHARACTER(LEN=MAX_NAME_LEN) :: SolverName,CostFile,UsedDataFile
+  CHARACTER(LEN=MAX_NAME_LEN) :: SolverName="Adjoint Solver"
+  CHARACTER(LEN=MAX_NAME_LEN) :: CostFile,UsedDataFile
   CHARACTER(LEN=MAX_NAME_LEN) :: CostSolName,VariableName,ObsFileName
   TYPE(Mesh_t),POINTER :: Mesh
   TYPE(Element_t),POINTER ::  Element
@@ -116,6 +117,13 @@ SUBROUTINE AdjointSSA_CostDiscSolver( Model,Solver,dt,TransientSimulation )
   save SolverName,CostSolName,CostFile,VariableName
   save ElementNodes
 
+   CALL Info(SolverName,'***********************',level=0)
+   CALL Info(SolverName,' This solver has been replaced by:',level=0)
+   CALL Info(SolverName,'   Adjoint_CostDiscSolver  ',level=0)
+   CALL Info(SolverName,' See documentation under:   ',level=0)
+   CALL Info(SolverName,'   elmerice/Solvers/Documentation   ',level=0)
+   CALL Info(SolverName,'***********************',level=0)
+   CALL FATAL(SolverName,' Use new solver !!')
 
   SolverParams => GetSolverParams()
   DIM=GetInteger(SolverParams ,'Problem Dimension',Found)

elmerice/Solvers/AdjointSSA/AdjointSSA_CostRegSolver.F90

@@ -111,6 +111,14 @@ SUBROUTINE AdjointSSA_CostRegSolver( Model,Solver,dt,TransientSimulation )
 
    WRITE(SolverName, '(A)') 'CostSolver_Regular'
 
+   CALL Info(SolverName,'***********************',level=0)
+   CALL Info(SolverName,' This solver has been replaced by:',level=0)
+   CALL Info(SolverName,'   Adjoint_CostRegSolver  ',level=0)
+   CALL Info(SolverName,' See documentation under:   ',level=0)
+   CALL Info(SolverName,'   elmerice/Solvers/Documentation   ',level=0)
+   CALL Info(SolverName,'***********************',level=0)
+   CALL FATAL(SolverName,' Use new solver !!')
+
   SolverParams => GetSolverParams()
 
 !! Dimension of the pb; ie with SSA we can be 1D or 2D on a 2D mesh, or 2D on a 3D mesh

elmerice/Solvers/AdjointSolver.F90

@@ -85,11 +85,20 @@ SUBROUTINE AdjointSolver( Model,Solver,dt,TransientSimulation )
   Real(KIND=dp) :: Unorm
   REAL(KIND=dp), allocatable :: FullMat(:,:)
   character(len=50) :: fo1
-  character(LEN=MAX_NAME_LEN) :: SolName,SolverName
+  character(LEN=MAX_NAME_LEN) :: SolName
+  character(LEN=MAX_NAME_LEN) :: SolverName="Adjoint Solver"
 
 
   save SolverName,Firsttime,SolverInd,STIFF,FORCE,ExtPressure,LoadVector,Alpha,Beta,SlipCoeff,w
 
+   CALL Info(SolverName,'***********************',level=0)
+   CALL Info(SolverName,' This solver has been replaced by:',level=0)
+   CALL Info(SolverName,'   Adjoint_LinearSolver  ',level=0)
+   CALL Info(SolverName,' See documentation under:   ',level=0)
+   CALL Info(SolverName,'   elmerice/Solvers/Documentation   ',level=0)
+   CALL Info(SolverName,'***********************',level=0)
+   CALL FATAL(SolverName,' Use new solver !!')
+
   DIM = CoordinateSystemDimension()
 
   StiffMatrix => Solver % Matrix
@@ -104,7 +113,6 @@ SUBROUTINE AdjointSolver( Model,Solver,dt,TransientSimulation )
 
   if (Firsttime) then
           Firsttime=.False.
-          SolverName = "Adjoint Solver"
           N = Solver % Mesh % MaxElementDOFs
           allocate(FORCE(  2*NSDOFs*N ),STIFF( 2*NSDOFs*N,2*NSDOFs*N ),ExtPressure(N), & 
                     SlipCoeff(3,N),LoadVector(4,N),Alpha(N),Beta(N),w(N))

elmerice/Solvers/CostSolver_Adjoint.F90

@@ -114,6 +114,15 @@ SUBROUTINE CostSolver_Adjoint( Model,Solver,dt,TransientSimulation )
 
     WRITE(SolverName, '(A)') 'CostSolver_Adjoint'
 
+   CALL Info(SolverName,'***********************',level=0)
+   CALL Info(SolverName,' This solver has been replaced by:',level=0)
+   CALL Info(SolverName,'   Adjoint_CostContSolver  ',level=0)
+   CALL Info(SolverName,'   Adjoint_CostRegSolver  ',level=0)
+   CALL Info(SolverName,' See documentation under:   ',level=0)
+   CALL Info(SolverName,'   elmerice/Solvers/Documentation   ',level=0)
+   CALL Info(SolverName,'***********************',level=0)
+   CALL FATAL(SolverName,' Use new solver !!')
+
 !!!!!!! Check for parallel run 
     Parallel = (ParEnv % PEs > 1)