Commit 8913fcd4 authored by Benoit Urruty's avatar Benoit Urruty

mesh_24

parent 09476e8a
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.interpolate import griddata
from matplotlib.colors import SymLogNorm
from matplotlib import cm
import geopandas as gpd
import seaborn as sns
from matplotlib.colors import LogNorm
from scipy import ndimage
from mpl_toolkits.axes_grid1 import make_axes_locatable
import glob
size=40
params = {'legend.fontsize': 'xx-large',
'figure.titlesize':'xx-large',
'figure.figsize': (20,8),
'axes.labelsize': size,
'axes.titlesize': size,
'xtick.labelsize': size*0.75,
'ytick.labelsize': size*0.75,
'axes.titlepad': 25}
plt.rcParams.update(params)
def stat(x):
mean=np.nanmean(x)
median=np.nanquantile(x,0.5)
quartile1=np.nanquantile(x,0.25)
quartile3=np.nanquantile(x,0.75)
std=np.nanstd(x)
return mean,median,quartile1,quartile3,std
def figure(x,y,z,sf,titlevariable,colorbartitle,param):
stati=stat(z)
fig=plt.figure(figsize=[20,12])
plt.rcParams.update({'font.size': 15})
ax=plt.gca()
sf.plot(ax=ax, color='k')
im=plt.pcolormesh(x,y,z,norm=SymLogNorm(linthresh=param[0], linscale=param[1] , vmin=param[2], vmax=param[3]),cmap=cm.seismic)
cbar=plt.colorbar(im)
plt.xlim(x.min(),x.max())
plt.ylim(y.min(),y.max())
cbar.set_label(colorbartitle)
plt.title(titlevariable + 'Bedmap2-bedmachine \n mean=' + str(stati[0]) + ' median=' +str(stati[1]) + ' quartile 25%=' +str(stati[2]) +'\n quartile 75%=' +str(stati[3]) + ' std=' +str(stati[4]))
#fig.savefig(titlevariable + '.png')
plt.show()
#plt.close('all')
#fig1=plt.figure(figsize=[20,12])
#plt.rcParams.update({'font.size': 15})
#sns.boxplot(z)
#plt.close('all')
def line_flow(flowline,x,y,z_bed1,z_bed2):
#interpolate from line
x_fl,y_fl=np.loadtxt(flowline + '.csv',skiprows=1,delimiter=',',usecols=(22,23),unpack=True)
x_temp=(x_fl-x.min())/1000
y_temp=(y_fl-y.min())/1000
#bedmachine resolution 500m
x_tmp=(x_fl-x.min())/500
y_tmp=(y_fl-y.min())/500
# profil along line interpolated
f_bed1 = ndimage.map_coordinates(np.array(z_bed1),[y_temp,x_temp],mode='nearest', order=1)
f_bed2 = ndimage.map_coordinates(np.array(z_bed2),[y_tmp,x_tmp],mode='nearest', order=1)
dist=np.zeros(len(x_fl))
for i in range(0,len(x_fl)-1):
dist[i+1]=dist[i]+np.sqrt((x_fl[i]-x_fl[i+1])**2+(y_fl[i]-y_fl[i+1])**2)/1000
return x_fl,y_fl,f_bed1,f_bed2,dist
# DISCLAIMER: This function is copied from https://github.com/nwhitehead/swmixer/blob/master/swmixer.py,
# which was released under LGPL.
def resample_by_interpolation(signal, input_fs, output_fs):
scale = output_fs / input_fs
# calculate new length of sample
n = round(len(signal) * scale)
# use linear interpolation
# endpoint keyword means than linspace doesn't go all the way to 1.0
# If it did, there are some off-by-one errors
# e.g. scale=2.0, [1,2,3] should go to [1,1.5,2,2.5,3,3]
# but with endpoint=True, we get [1,1.4,1.8,2.2,2.6,3]
# Both are OK, but since resampling will often involve
# exact ratios (i.e. for 44100 to 22050 or vice versa)
# using endpoint=False gets less noise in the resampled sound
resampled_signal = np.interp(
np.linspace(0.0, 1.0, n, endpoint=False), # where to interpret
np.linspace(0.0, 1.0, len(signal), endpoint=False), # known positions
signal, # known data points
)
return resampled_signal
# chargement des données
sf = gpd.read_file(
"/home/urrutyb/Documents/PhD_TiPACCS/topo_antarctica/InSAR_GL_Antarctica_v02/InSAR_GL_Antarctica_v02.shp")
sf = sf.to_crs({'init': 'epsg:3031'})
bedmachine = xr.load_dataset(
"/home/urrutyb/Documents/PhD_TiPACCS/topo_antarctica/BedMachineAntarctica_2019-11-05_v01.nc")
figure='/home/urrutyb/Documents/PhD_TiPACCS/topo_antarctica/'
# selection des données
index = np.arange(0, len(bedmachine.x), 2)
bedmachinebed = np.flip(bedmachine.bed[index, index], 0)
x = np.array(bedmachinebed.x)
y = np.array(bedmachinebed.y)
del bedmachine,bedmachinebed
mask_out_error=np.loadtxt('../../../backup.txt',delimiter=',')
# mask_out_error[index]=np.nan
# mask_out_error[index2]=np.nan
# np.savetxt('backup.txt',mask_out_error,delimiter=',')
fig=plt.figure(figsize=[12,10])
ax=plt.gca()
sf.plot(ax=ax, color='k')
im=plt.pcolormesh(x,y,mask_out_error)
cbar=plt.colorbar(im)
cbar.set_label('1 : overlap - 0 : outside')
plt.axis('off')
# plt.title('link between the two beds +/- uncertainty')
plt.tight_layout()
plt.savefig(figure + 'uncertainty.png',dpi=300)
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.interpolate import griddata
from matplotlib.colors import SymLogNorm
from matplotlib import cm
import geopandas as gpd
import seaborn as sns
from matplotlib.colors import LogNorm
from scipy import ndimage
from mpl_toolkits.axes_grid1 import make_axes_locatable
import glob
size=40
params = {'legend.fontsize': 'xx-large',
'figure.titlesize':'xx-large',
'figure.figsize': (20,8),
'axes.labelsize': size,
'axes.titlesize': size,
'xtick.labelsize': size*0.5,
'ytick.labelsize': size*0.5,
'axes.titlepad': 25}
plt.rcParams.update(params)
def stat(x):
mean=np.nanmean(x)
median=np.nanquantile(x,0.5)
quartile1=np.nanquantile(x,0.25)
quartile3=np.nanquantile(x,0.75)
std=np.nanstd(x)
return mean,median,quartile1,quartile3,std
def figure(x,y,z,sf,titlevariable,colorbartitle,param):
stati=stat(z)
fig=plt.figure(figsize=[20,12])
plt.rcParams.update({'font.size': 15})
ax=plt.gca()
sf.plot(ax=ax, color='k')
im=plt.pcolormesh(x,y,z,norm=SymLogNorm(linthresh=param[0], linscale=param[1] , vmin=param[2], vmax=param[3]),cmap=cm.seismic)
cbar=plt.colorbar(im)
plt.xlim(x.min(),x.max())
plt.ylim(y.min(),y.max())
cbar.set_label(colorbartitle)
plt.title(titlevariable + 'Bedmap2-bedmachine \n mean=' + str(stati[0]) + ' median=' +str(stati[1]) + ' quartile 25%=' +str(stati[2]) +'\n quartile 75%=' +str(stati[3]) + ' std=' +str(stati[4]))
#fig.savefig(titlevariable + '.png')
plt.show()
#plt.close('all')
#fig1=plt.figure(figsize=[20,12])
#plt.rcParams.update({'font.size': 15})
#sns.boxplot(z)
#plt.close('all')
def line_flow(flowline,x,y,z_bed1,z_bed2):
#interpolate from line
x_fl,y_fl=np.loadtxt(flowline + '.csv',skiprows=1,delimiter=',',usecols=(22,23),unpack=True)
x_temp=(x_fl-x.min())/1000
y_temp=(y_fl-y.min())/1000
#bedmachine resolution 500m
x_tmp=(x_fl-x.min())/500
y_tmp=(y_fl-y.min())/500
# profil along line interpolated
f_bed1 = ndimage.map_coordinates(np.array(z_bed1),[y_temp,x_temp],mode='nearest', order=1)
f_bed2 = ndimage.map_coordinates(np.array(z_bed2),[y_tmp,x_tmp],mode='nearest', order=1)
dist=np.zeros(len(x_fl))
for i in range(0,len(x_fl)-1):
dist[i+1]=dist[i]+np.sqrt((x_fl[i]-x_fl[i+1])**2+(y_fl[i]-y_fl[i+1])**2)/1000
return x_fl,y_fl,f_bed1,f_bed2,dist
# DISCLAIMER: This function is copied from https://github.com/nwhitehead/swmixer/blob/master/swmixer.py,
# which was released under LGPL.
def resample_by_interpolation(signal, input_fs, output_fs):
scale = output_fs / input_fs
# calculate new length of sample
n = round(len(signal) * scale)
# use linear interpolation
# endpoint keyword means than linspace doesn't go all the way to 1.0
# If it did, there are some off-by-one errors
# e.g. scale=2.0, [1,2,3] should go to [1,1.5,2,2.5,3,3]
# but with endpoint=True, we get [1,1.4,1.8,2.2,2.6,3]
# Both are OK, but since resampling will often involve
# exact ratios (i.e. for 44100 to 22050 or vice versa)
# using endpoint=False gets less noise in the resampled sound
resampled_signal = np.interp(
np.linspace(0.0, 1.0, n, endpoint=False), # where to interpret
np.linspace(0.0, 1.0, len(signal), endpoint=False), # known positions
signal, # known data points
)
return resampled_signal
# chargement des données
sf = gpd.read_file(
"/home/urrutyb/Documents/PhD_TiPACCS/topo_antarctica/InSAR_GL_Antarctica_v02/InSAR_GL_Antarctica_v02.shp")
sf = sf.to_crs({'init': 'epsg:3031'})
bedmachine = xr.load_dataset(
"/home/urrutyb/Documents/PhD_TiPACCS/topo_antarctica/BedMachineAntarctica_2019-11-05_v01.nc")
figure='/home/urrutyb/Documents/PhD_TiPACCS/topo_antarctica/'
# selection des données
index = np.arange(0, len(bedmachine.x), 2)
bedmachinebed = np.flip(bedmachine.bed[index, index], 0)
x = np.array(bedmachinebed.x)
y = np.array(bedmachinebed.y)
del bedmachine,bedmachinebed
mask_out_error=np.loadtxt('../../../backup.txt',delimiter=',')
# mask_out_error[index]=np.nan
# mask_out_error[index2]=np.nan
# np.savetxt('backup.txt',mask_out_error,delimiter=',')
fig=plt.figure(figsize=[24,20])
ax=plt.gca()
sf.plot(ax=ax, color='k')
im=plt.pcolormesh(x,y,mask_out_error)
cbar=plt.colorbar(im)
cbar.set_label('1 : overlap - 0 : outside')
# plt.title('link between the two beds +/- uncertainty')
plt.tight_layout()
plt.savefig(figure + 'uncertainty.png',dpi=300)
s# gmsh command line
gmsh -1 -2 -format msh2 Contour.geo
# command to generate mesh for elmerice
ElmerGrid 14 2 Contour.msh -autoclean
# command to generate vtu from msh file
ElmerGrid 14 5 Contour.msh -autoclean
This diff is collapsed.
This source diff could not be displayed because it is too large. You can view the blob instead.
This diff is collapsed.
This diff is collapsed.
!##################################
! 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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! .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 = "ScatteredInter"
Procedure = "Scattered2DDataInterpolator" "Scattered2DDataInterpolator"
Variable = -nooutput dumy
Bounding Box dx = Real 20000.0
!!!!! NNI or linear (nn-c library)
Variable 1 = File "thickness"
Target Variable 1 = String "H"
Variable 1 data file = File "$TOPOGRAPHY_DATA$"
Variable 1 method = String "li"
Variable 2 = File "norm"
Target Variable 2 = String "Vobs"
Variable 2 data file = File "$VELOCITY_DATA$"
Exported Variable 1 = H
Exported Variable 2 = "Vobs"
End
Solver 2
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 3
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 4
Equation = "Metric1"
Procedure = "ElmerIce_MeshAdapt2D" "MMG2D_MetricAniso"
Metric Variable Name = String "Mu"
Hessian Variable Name = String "ddx1"