#Test Name: EarthSlr_rotationalFeedback
import numpy as np

from socket import gethostname

from gmshplanet import *
from gmtmask import *
from lovenumbers import *
from model import *
from paterson import *
from solve import *


#mesh earth:
md = model()
md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 1000.)  #1000 km resolution mesh

#parameterize solidearth solution:
#solidearth loading:  {{{
md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, ))
md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, ))
md.dsl.global_average_thermosteric_sea_level_change=np.zeros((2, ))
md.dsl.sea_surface_height_change_above_geoid=np.zeros((md.mesh.numberofvertices+1, ))
md.dsl.sea_water_pressure_change_at_sea_floor=np.zeros((md.mesh.numberofvertices+1, ))

#antarctica
#Access every element in lat using the indices in elements
# - 1 to convert to base 0 indexing, 1 (not 2, in matlab) to sum over rows
late = md.mesh.lat[md.mesh.elements - 1].sum(axis=1)/ 3
longe = md.mesh.long[md.mesh.elements - 1].sum(axis=1)/ 3
pos = np.intersect1d(np.array(np.where(late < -75)[0]), np.array(np.where(longe < 0)[0]))
md.solidearth.surfaceload.icethicknesschange[pos] = -1

#elastic loading from love numbers:
md.solidearth.lovenumbers = lovenumbers('maxdeg', 1000)
#}}}

#mask:  {{{
mask = gmtmask(md.mesh.lat, md.mesh.long)
icemask = np.ones(md.mesh.numberofvertices)
pos = np.where(mask == 0)[0]
icemask[pos] = -1
pos = np.where(mask[md.mesh.elements - 1].sum(axis=1) < 3)[0]
icemask[md.mesh.elements[pos, :] - 1] = -1
md.mask.ice_levelset = icemask
md.mask.ocean_levelset = -icemask

#make sure that the elements that have loads are fully grounded
pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0]
md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1

#make sure wherever there is an ice load, that the mask is set to ice:
#pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] # Do we need to do this twice?
md.mask.ice_levelset[md.mesh.elements[pos, :] - 1] = -1
# }}}

# use model representation of ocen area (not the true area)
md.solidearth.settings.ocean_area_scaling = 0

#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, ))

#Miscellaneous
md.miscellaneous.name = 'test2003'

#Solution parameters
md.solidearth.settings.reltol = np.nan
md.solidearth.settings.abstol = 1e-3
md.solidearth.settings.computesealevelchange = 1

#eustatic + rigid + elastic run:
md.solidearth.settings.rigid = 1
md.solidearth.settings.elastic = 1
md.solidearth.settings.rotation = 0
md.cluster = generic('name', gethostname(), 'np', 3)
#md.verbose = verbose('111111111')
md = solve(md, 'Sealevelrise')
SnoRotation = md.results.SealevelriseSolution.Sealevel

#eustatic + rigid + elastic + rotation run:
md.solidearth.settings.rigid = 1
md.solidearth.settings.elastic = 1
md.solidearth.settings.rotation = 1
md.cluster = generic('name', gethostname(), 'np', 3)
#md.verbose = verbose('111111111')
md = solve(md, 'Sealevelrise')
SRotation = md.results.SealevelriseSolution.Sealevel

#Fields and tolerances to track changes
field_names = ['noRotation', 'Rotation']
field_tolerances = [1e-13, 1e-13]
field_values = [SnoRotation, SRotation]
