source: issm/trunk/test/Par/Pig.py@ 16560

import os.path
import inspect
import netCDF4
import numpy
from verbose import verbose
from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d
from paterson import paterson 
from SetMarineIceSheetBC import SetMarineIceSheetBC

#Start defining model parameters here

#Geometry and observation
f = netCDF4.Dataset('../Data/Pig.nc','r')
x         = numpy.reshape(f.variables['x'][:],(-1))
y         = numpy.reshape(f.variables['y'][:],(-1))
vx_obs    = f.variables['vx_obs'][:]
vy_obs    = f.variables['vy_obs'][:]
index     = f.variables['index'][:]
surface   = f.variables['surface'][:]
thickness = f.variables['thickness'][:]
f.close()

[md.inversion.vx_obs]   =InterpFromMeshToMesh2d(index,x,y,vx_obs,md.mesh.x,md.mesh.y)
[md.inversion.vy_obs]   =InterpFromMeshToMesh2d(index,x,y,vy_obs,md.mesh.x,md.mesh.y)
[md.geometry.surface]  =InterpFromMeshToMesh2d(index,x,y,surface,md.mesh.x,md.mesh.y)
[md.geometry.thickness]=InterpFromMeshToMesh2d(index,x,y,thickness,md.mesh.x,md.mesh.y)
md.geometry.bed=md.geometry.surface-md.geometry.thickness
md.initialization.vx=md.inversion.vx_obs
md.initialization.vy=md.inversion.vy_obs
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))
md.initialization.temperature=md.initialization.temperature

#Damage
md.damage.D=numpy.zeros((md.mesh.numberofvertices,1))

#Friction
md.friction.coefficient=50.*numpy.ones((md.mesh.numberofvertices,1))
md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
md.friction.p=numpy.ones((md.mesh.numberofelements,1))
md.friction.q=numpy.ones((md.mesh.numberofelements,1))

#Numerical parameters
md.stressbalance.viscosity_overshoot=0.3
md.masstransport.stabilization=1.
md.verbose=verbose(0)
md.settings.waitonlock=30
md.timestepping.time_step=1.
md.timestepping.final_time=2.
md.stressbalance.restol=0.05
md.stressbalance.reltol=1.
md.steadystate.reltol=1.
md.stressbalance.abstol=float('nan')

#Boundary conditions:
md=SetMarineIceSheetBC(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]