source: issm/trunk-jpl/test/NightlyRun/test2113.py@ 24261

#Test Name: Esa2Dsurface
#Northern hemisphere example for north - south, east - west components of horiz motion
#Same as test2111.m except that AIS is assumed to have located in Northern Hemisphere

import numpy as np
from model import *
from socket import gethostname
from solve import *
from roundmesh import *
from love_numbers import *
from paterson import *

#mesh ais: {{{
md = model()
md = triangle(md, '../Exp/Ais.exp', 200000)  # max element size
# }}}
#define load: {{{
md.esa.deltathickness = np.zeros((md.mesh.numberofelements, ))
disc_radius = 500  # km
index = md.mesh.elements
x_element = np.mean(md.mesh.x[index - 1], 1) - 1.0e6
y_element = np.mean(md.mesh.y[index - 1], 1) + 1.0e6
rad_dist = np.sqrt(x_element**2 + y_element**2) / 1000  # radial distance in km
md.esa.deltathickness[np.where(rad_dist <= disc_radius)] = -1  # 1 m water withdrawl
# }}}
#love numbers: {{{
nlov = 10000  # horizontal displacements do not work for low degree truncation, e.g., 101
md.esa.love_h = np.array(love_numbers('h', 'CF'))
md.esa.love_h = np.resize(md.esa.love_h, nlov + 1)
md.esa.love_l = np.array(love_numbers('l', 'CF'))
md.esa.love_l = np.resize(md.esa.love_l, nlov + 1)
# }}}
#mask:  {{{
#make sure wherever there is an ice load, that the mask is set to ice:
md.mask.ice_levelset = np.ones((md.mesh.numberofvertices, ))
pos = np.where(md.esa.deltathickness)
md.mask.ice_levelset[md.mesh.elements[pos, :]] = -1

#is ice grounded?
md.mask.groundedice_levelset = -np.ones((md.mesh.numberofvertices, ))
pos = np.where(md.mask.ice_levelset <= 0)
md.mask.groundedice_levelset[pos] = 1
# }}}
#geometry:  {{{
di = md.materials.rho_ice / md.materials.rho_water
md.geometry.thickness = np.ones((md.mesh.numberofvertices, ))
md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices, ))
md.geometry.base = md.geometry.surface - md.geometry.thickness
md.geometry.bed = md.geometry.base
# }}}
#materials:  {{{
md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices, ))
md.materials.rheology_B = paterson(md.initialization.temperature)
md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, ))
# }}}
#additional parameters, miscellaneous: {{{
md.miscellaneous.name = 'test2113'
md.esa.degacc = 0.01
md.esa.hemisphere = 1  # AIS is placed in Northern Hemisphere
# }}}

#solve esa: {{{
md.esa.requested_outputs = ['EsaUmotion', 'EsaNmotion', 'EsaEmotion']
md.cluster = generic('name', gethostname(), 'np', 3)
md.verbose = verbose('111111111')
md = solve(md, 'Esa')
# }}}
#Fields and tolerances to track changes: {{{
field_names = ['EsaUmotion', 'EsaNmotion', 'EsaEmotion']
field_tolerances = [1e-13, 1e-13, 1e-13]
field_values = [md.results.EsaSolution.EsaUmotion,
                md.results.EsaSolution.EsaNmotion,
                md.results.EsaSolution.EsaEmotion]
# }}}