Commit a7ccd88f authored by Benoit Urruty's avatar Benoit Urruty

inversion

parent 95386375
UTF-8
\ No newline at end of file
PROJCS["WGS_84_Antarctic_Polar_Stereographic",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Stereographic_South_Pole"],PARAMETER["standard_parallel_1",-71],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]
\ No newline at end of file
PROJCS["WGS 84 / Antarctic Polar Stereographic",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",-71],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3031"]]
PROJCS["WGS_84_Antarctic_Polar_Stereographic",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Stereographic_South_Pole"],PARAMETER["standard_parallel_1",-71],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]
\ No newline at end of file
PROJCS["WGS 84 / Antarctic Polar Stereographic",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",-71],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3031"]]
UTF-8
\ No newline at end of file
PROJCS["WGS_84_Antarctic_Polar_Stereographic",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Stereographic_South_Pole"],PARAMETER["standard_parallel_1",-71],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]
\ No newline at end of file
PROJCS["WGS 84 / Antarctic Polar Stereographic",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",-71],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3031"]]
# Inversion of the viscosity and the friction parameter for the Ronne-Fichner Ice Shelf
## Creation of the mesh
## Initialization
We need to interpolate the data on the mesh. These data contains :
- The bed topography, the thickness (both from BedMachine)
- The surface mass balance (SMB from a the MAR model)
- The mean viscosity which correspond to a first guess on $\mu$ (from Paterson)
- The surface velocity (from MEaSUREs)
The models is now able to compute others variables.
- The surface topography
- The groundedmask
- The basal mass balance based on the Pollard & De Conto parametrization
- The nodal gradient
- The first guess on $\alpha$ (the friction parameter)
## Inversion
The goal of the inversion is to find the viscosity and the friction parameter which will permit
to approach as best as possible the observed surface velocity.
### SSA
The SSA is compute with the first guess on the value to optimize. We
solve the equation :
$$
\frac{\partial}{\partial x}(2\eta h (2 \frac{\partial u}{\partial x} + \eta \frac{\partial v}{\partial y})) + \frac{\partial}{\partial y}(\eta h (\frac{\partial u}{\partial y} + \eta \frac{\partial v}{\partial x})) - \beta^2u =\rho_igh\frac{\partial s}{\partial x}
$$
$$
\frac{\partial}{\partial y}(2\eta h (2 \frac{\partial v}{\partial y} + \eta \frac{\partial u}{\partial x})) + \frac{\partial}{\partial x}(\eta h (\frac{\partial v}{\partial x} + \eta \frac{\partial u}{\partial y})) - \beta^2v =\rho_igh\frac{\partial s}{\partial y}
$$
The effective viscosity :
$$
\eta = \frac{1}{2}A^{-1/n_{\epsilon_e}(1-n)/n}
$$
The value of the velocity u and v are exported to be used in the cost function.
### The cost function
The cost function is the comparison of the compute surface velocity
and the observed one.
$$
J = \int{\frac{1}{2}((u_{SSA} - u_{Obs})^2+(v_{SSA} - v_{Obs})^2)} d\Omega
$$
Derivative of the cost function
$$
velocityb = \int{((u_{SSA} - u_{Obs})\partial u+(v_{SSA} - v_{Obs})\partial v)} d\Omega
$$
### The adjoint method
The adjoint problem is a numerical method to compute the gradient in a numerical optimization problem. The adjoint for a linear system **A**x = F is define as :
$$
\boldsymbol{A}^T b = xb
$$
### The gradient
We need to compute the gradient for each parameters we have to optimize. Theoretically, we have the right solution when the gradient is equal to zero but we never obtains this value. So we define a lower limit which is define a the gradient is enough low to have a good result.
$$
\nabla J= \left\{ \begin{array}{ll}
\frac{\partial J}{\partial Beta}=\\
\frac{\partial J}{\partial Eta}=
\end{array}
\right\}
$$
### The regularization terms
The regularization can be done one the first derivative. We gave a weight $\lambda$ to this the value computed. This weight have to be choice according the
$$
J_{reg} = \int_{\Omega} 0.5 (|dV/dx|)^2 d\Omega
\\with\ V \ the\ nodal\ variable
$$
$$
J_{reg} = \int_{\Omega} 0.5 ((V-V^{prior})/s^2)^2 d\Omega
\\with\ V \ the\ nodal\ variable, V^{prior}\ is \ the\ estimate\ and\ s^2\ the\ variance.
$$
### The optimization M1QN3
The optimization will take the gradient and try to found the best choice by doing simulations to find the best direction and will pass to the next iteration.
\ No newline at end of file
// This a a geo file created using the python script Makegeo.py // // This a a geo file created using the python script Makegeo.py //
Mesh.Algorithm=5; Mesh.Algorithm=5;
// To controle the element size, one can directly modify the lc value in the geo file // // To controle the element size, one can directly modify the lc value in the geo file //
lc = 1000 ; lc = 10000 ;
Point(1) = { -1461250.0, 782211.997539443, 0.0, lc}; Point(1) = { -1461250.0, 782211.997539443, 0.0, lc};
Point(2) = { -1461250.0, 782750.0, 0.0, lc}; Point(2) = { -1461250.0, 782750.0, 0.0, lc};
Point(3) = { -1460750.0, 782750.0, 0.0, lc}; Point(3) = { -1460750.0, 782750.0, 0.0, lc};
...@@ -2623,3 +2623,12 @@ Plane Surface(1310) = {1309}; ...@@ -2623,3 +2623,12 @@ Plane Surface(1310) = {1309};
Physical Curve(1311) = {1:973}; Physical Curve(1311) = {1:973};
Physical Curve(1312) = {973:1308}; Physical Curve(1312) = {973:1308};
Physical Surface(1313) = {1310}; Physical Surface(1313) = {1310};
Field[1] = MathEval;
Field[1].F = Sprintf("5000.0");
Background Field = 1;
Mesh.CharacteristicLengthFromCurvature = 0;
Mesh.CharacteristicLengthExtendFromBoundary = 0;
Mesh.CharacteristicLengthFromPoints = 1;
This diff is collapsed.
...@@ -18,7 +18,7 @@ ...@@ -18,7 +18,7 @@
{ {
"data": { "data": {
"text/plain": [ "text/plain": [
"<matplotlib.collections.PathCollection at 0x7f2c3f128bd0>" "<matplotlib.collections.PathCollection at 0x7fb4df93cb10>"
] ]
}, },
"execution_count": 1, "execution_count": 1,
...@@ -67,7 +67,7 @@ ...@@ -67,7 +67,7 @@
{ {
"data": { "data": {
"text/plain": [ "text/plain": [
"<matplotlib.collections.PathCollection at 0x7f2c3c8bac50>" "<matplotlib.collections.PathCollection at 0x7fb4df93c950>"
] ]
}, },
"execution_count": 2, "execution_count": 2,
...@@ -117,7 +117,7 @@ ...@@ -117,7 +117,7 @@
{ {
"data": { "data": {
"text/plain": [ "text/plain": [
"[<matplotlib.lines.Line2D at 0x7f2c3c7e4d90>]" "[<matplotlib.lines.Line2D at 0x7fb4df7e87d0>]"
] ]
}, },
"execution_count": 4, "execution_count": 4,
...@@ -145,7 +145,27 @@ ...@@ -145,7 +145,27 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 21, "execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-194344.40573864247"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mesh[:,1].min()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
......
!##################################
! INCLUDE file for MESH_OPTIM.sif
!###################################
!name of the RUN : optimised mesh will be saved in MESH_$RUN
$RUN="1"
! the name of the mesh to optimize
$MESH_IN="MESH"
!##########################################################
!##########################################################
!## Mesh criteria
!##########################################################
!##########################################################
$IMIN=3
$IMAX=10
$Tol=0.01
!##########################################################
!##########################################################
!## Mesh criteria
!##########################################################
!##########################################################
! Tolerated errors on U and H
$U_err=15.0
$H_err=20.0
! Minimal and maximal mesh size
$Hmin=500.0
!Maximal mesh size is function to the disance to the margin
$MarginExtent=40.0e3
! Maximal mesh size outside the margin
$HMaxOUT=50.0e3
! Maximal mesh size within the margin
$HMaxIN=4.0e3
!##########################################################
!##########################################################
!###### DATA FILES
!##########################################################
!##########################################################
$TOPOGRAPHY_DATA="/home/urrutyb/Documents/PhD_TiPACCS/topo_antarctica/BedMachineAntarctica_2019-11-05_v01_m.nc"
$VELOCITY_DATA="/home/urrutyb/Documents/PhD_TiPACCS/topo_antarctica/ronne_bassin_velocity_norm.nc"
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! .sif file for greenland mesh adaptation
! mesh is adapted to equidistribute the interpolation error
! of the thickness and velocity observation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! FOR DEFAULT USE UPDATE PARAMETERS IN MESH_OPTIM.IN
include MESH_OPTIM.IN
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Header
Mesh DB "." "$MESH_IN$"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Constants
Hmax margin extent = Real $MarginExtent
Hmax within margin = Real $HMaxIN
Hmax outside margin = Real $HMaxOUT
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Simulation
Coordinate System = Cartesian 2D
Simulation Type = Steady
Steady State Min Iterations = $IMAX
Steady State Max Iterations = $IMAX
max output level = 3
Post File = "OPTIM_$RUN$.vtu"
! OutPut File = "OPTIM_$RUN$.result"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body 1
Equation = 1
Body Force = 1
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body Force 1
Mu Hmin = Real $Hmin
Mu Hmax = Variable Distance
REAL procedure "USF_PCond" "Hmax"
Mu err = Real $U_err
Mh Hmin = Real $Hmin
Mh Hmax = Variable Distance
REAL procedure "USF_PCond" "Hmax"
Mh err = Real $H_err
Distance = Real 0.0
Distance Condition = Variable H
Real procedure "USF_PCond" "DistanceCond"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Solver 1
Equation = Reader
Procedure = "GridDataReader" "GridDataReader"
Variable = -nooutput dumy
!---- NOTE: File is case sensitive, String is not!
Filename = File "$TOPOGRAPHY_DATA$"
X Dim Name = String "x"
Y Dim Name = String "y"
!---
Variable 1 = File "thickness"
Target Variable 1 = String "H"
Exported Variable 1 = H
End
Solver 2
Equation = Reader2
Procedure = "GridDataReader" "GridDataReader"
Variable = -nooutput dumy
!---- NOTE: File is case sensitive, String is not!
Filename = File "$VELOCITY_DATA$"
X Dim Name = String "x"
Y Dim Name = String "y"
!--- Interpolation variables
Variable 1 = File "norm"
Target Variable 1 = String "Vobs"
Exported Variable 1 = "Vobs"
End
Solver 3
Equation = "Distance"
Variable = distance
Procedure = "DistanceSolve" "DistanceSolver1"
Optimize Bandwidth = logical false
End
! Compute the metric associated with f1
! Compute the metric associated with f1
! 1- compute projected gradient of f1
Solver 4
Equation = "Nodal Gradient 1"
Variable = -nooutput "Gradient1"
Variable DOFs = 2
Procedure = "ElmerIce_MeshAdapt2D" "Compute2DNodalGradient"
Variable Name = string "vobs"
End
! 2- compute:
! - the hessian matrix by solving a diffusion equation:
! - the metric tensor
Solver 5
Equation = "Metric1"
Procedure = "ElmerIce_MeshAdapt2D" "MMG2D_MetricAniso"
Metric Variable Name = String "Mu"
Hessian Variable Name = String "ddx1"
Gradient Name = String "Gradient1"
Diffusivity = Real 0.5
Linear System Solver = Direct
Linear System Direct Method = umfpack
Linear System Refactorize = False
End
! Compute the metric associated with f2
Solver 6
Equation = "Nodal Gradient 2"
Variable = -nooutput "Gradient2"
Variable DOFs = 2
Procedure = "ElmerIce_MeshAdapt2D" "Compute2DNodalGradient"
Variable Name = string "h"
End
Solver 7
Equation = "Metric2"
Procedure = "ElmerIce_MeshAdapt2D" "MMG2D_MetricAniso"
Metric Variable Name = String "Mh"
Hessian Variable Name = String "ddx2"
Gradient Name = String "Gradient2"
Diffusivity = Real 0.5
Linear System Solver = Direct
Linear System Direct Method = umfpack
Linear System Refactorize = False
End
!! do the intersection of M1 and M2
Solver 8
Equation = "Metric"
Procedure = "ElmerIce_MeshAdapt2D" "MMG2D_MetricIntersect"
Metric Variable Name = String "M1M2"
Metric 1 Variable Name = String "Mh"
Metric 2 Variable Name = String "Mu"
End
Solver 9
Equation = SaveScalars
Procedure = "SaveData" "SaveScalars"
Filename = "f_$RUN$.dat"
Show Norm Index = Integer 2
Variable 1 = "Time"
Operator 2 = nodes
Variable 3 = "H"
Operator 3 = "int"
Variable 4 = "Vobs"
Operator 4 = "int"
End
!! Anisotropic mesh adaptation using the MMG library
Solver 10
!! mandatory else Model % Mesh % Changed reset to .FALSE. in coupled simulations
Exec Solver = after timestep
Equation = "MMG"
Procedure = "ElmerIce_MeshAdapt2D" "MMG2DSolver"
Output file name = "MESH_$RUN$"
Metric Variable Name = String "M1M2"
Angle detection = Real 5.0
Increment Mesh Number = logical false
Release previous mesh = Logical True
Steady State Convergence Tolerance = Real $Tol
Steady State Min Iterations = INTEGER $IMIN
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Equation 1
Active Solvers(10) = 1 2 3 4 5 6 7 8 9 10
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Boundary Condition 1
Target Boundaries(1) = 1
End
LIBS=USF_PCond
all: $(LIBS)
%: %.F90
elmerf90 $^ -o $@
clean:
rm -f $(LIBS)
This directory contains standard configuration files to adapt an initial mesh of the Greenland Ice Sheet:
see http://elmerice.elmerfem.org/wiki/doku.php?id=eis:greenland#present
USAGE:
- Get an initial mesh : e.g. http://elmerice.elmerfem.org/wiki/lib/exe/fetch.php?media=eis:greenland:present:greenland_mesh_v0.tar.gz
- Update default parameteris in MESH_OPTIM.IN
- (OPTIONAL - Adavnaced users): update default .sif file MESH_OPTIM.sif
- compile user-function: make -f Makefile.IN
- run the mesh optimisation: ElmerSolver
FUNCTION DistanceCond(Model,nodenumber,VarIn) RESULT(VarOut)
USE DefUtils
implicit none
!-----------------
TYPE(Model_t) :: Model
INTEGER :: nodenumber
REAL(kind=dp) :: VarIn,VarOut
IF (VarIn.LT.0.1) THEN
VarOut = +1.0
ELSE
VarOut = -1.0
END IF
End FUNCTION DistanceCond
FUNCTION HMax(Model,nodenumber,VarIn) RESULT(VarOut)
USE DefUtils
implicit none
!-----------------
TYPE(Model_t) :: Model
INTEGER :: nodenumber
REAL(kind=dp) :: VarIn !Distance
REAL(kind=dp) :: VarOut
REAL(kind=dp),SAVE :: Extent,HMaxIN,HMaxOut
LOGICAL,SAVE :: FirstVisit=.TRUE.
IF (FirstVisit) THEN
Extent=ListGetConstReal( Model % Constants,'Hmax margin extent',UnFoundFatal=.TRUE.)
HMaxIN=ListGetConstReal( Model % Constants,'Hmax within margin',UnFoundFatal=.TRUE.)
HMaxOut=ListGetConstReal( Model % Constants,'Hmax outside margin',UnFoundFatal=.TRUE.)
END IF
IF ((VarIn.GT.1.0).AND.(VarIn.LT.Extent)) THEN
VarOut = HMaxIN
ELSE
VarOut = HMaxOut
END IF
End FUNCTION HMax