source: issm/trunk-jpl/test/Par/SquareSheetConstrained.py@ 14826

import os.path
import netCDF4
import numpy
import inspect
from verbose import *
from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d
from paterson import *
from SetIceSheetBC import *

#Start defining model parameters here

#Geometry
hmin=300.
hmax=1000.
ymin=numpy.min(md.mesh.y)
ymax=numpy.max(md.mesh.y)

md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y.reshape(-1,1)-ymin)/(ymax-ymin)
md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness+20.
md.geometry.surface=md.geometry.bed+md.geometry.thickness

#Initial velocity 
f = netCDF4.Dataset('../Data/SquareSheetConstrained.nc','r')
x     = numpy.reshape(f.variables['x'][:],(-1))
y     = numpy.reshape(f.variables['y'][:],(-1))
vx    = f.variables['vx'][:]
vy    = f.variables['vy'][:]
index = f.variables['index'][:]
f.close()

[md.initialization.vx]=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y)
[md.initialization.vy]=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y)
md.initialization.vz=numpy.zeros((md.mesh.numberofvertices,1))
md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1))

#Materials
md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1))
md.materials.rheology_B=paterson(md.initialization.temperature)
md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))

#Friction
pos=numpy.nonzero(md.mask.elementonfloatingice)
md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices,1))
md.friction.coefficient[md.mesh.elements[pos,:].astype(int)-1]=0.
md.friction.p=numpy.ones((md.mesh.numberofelements,1))
md.friction.q=numpy.ones((md.mesh.numberofelements,1))

#Numerical parameters
md.diagnostic.viscosity_overshoot=0.0
md.prognostic.stabilization=1.
md.thermal.stabilization=1.
md.verbose=verbose(0)
md.settings.waitonlock=30
md.diagnostic.restol=0.05
md.steadystate.reltol=0.05
md.diagnostic.reltol=0.05
md.diagnostic.abstol=float('NaN')
md.timestepping.time_step=1.
md.timestepping.final_time=3.

#GIA:
md.gia.lithosphere_thickness=100.*numpy.ones((md.mesh.numberofvertices,1)); # in km
md.materials.lithosphere_shear_modulus=6.7*10**10;                          # in Pa
md.materials.lithosphere_density=3.32;                                      # in g/cm^-3
md.materials.mantle_shear_modulus=1.45*10**11;                              # in Pa
md.materials.mantle_viscosity=1.*10**21;                                    # in Pa.s
md.materials.mantle_density=3.34;                                           # in g/cm^-3

#Boundary conditions:
md=SetIceSheetBC(md)

#Change name so that no test have the same name
if len(inspect.stack()) > 2:
        md.miscellaneous.name = os.path.basename(inspect.stack()[2][1]).split('.')[0]