#Test Name: SquareNoDynHydrologyDCSmbCoupled
import numpy as np
from model import *
from socket import gethostname
from triangle import triangle
from setmask import setmask
from parameterize import parameterize
from solve import solve
from generic import generic
from SMBgradientscomponents import SMBgradientscomponents

md = triangle(model(), '../Exp/Square.exp', 100000.)
md = setmask(md, '', '')
md = parameterize(md, '../Par/SquareNoDyn.py')
md.cluster = generic('name', gethostname(), 'np', 1)

md.transient.ishydrology = True
md.transient.issmb = True
md.hydrology = hydrologydc()
md.hydrology = md.hydrology.initialize(md)
md.smb = SMBgradientscomponents()

md.hydrology.isefficientlayer = 1

md.hydrology.sedimentlimit_flag = 1
md.hydrology.sedimentlimit = 400.0
md.hydrology.sediment_transmitivity = 3.0 * np.ones((md.mesh.numberofvertices))
md.hydrology.mask_thawed_node = np.ones((md.mesh.numberofvertices))

md.hydrology.mask_eplactive_node = np.zeros((md.mesh.numberofvertices))
md.hydrology.epl_conductivity = 3.
md.hydrology.epl_initial_thickness = 20
md.hydrology.epl_colapse_thickness = 1.0e-3
md.hydrology.epl_thick_comp = 0
md.hydrology.epl_max_thickness = 1

md.hydrology.spcsediment_head = np.nan * np.ones((md.mesh.numberofvertices))
md.hydrology.spcepl_head = np.nan * np.ones((md.mesh.numberofvertices))

md.initialization.sediment_head = np.zeros((md.mesh.numberofvertices))
md.initialization.epl_head = np.zeros((md.mesh.numberofvertices))
md.initialization.epl_thickness = np.ones((md.mesh.numberofvertices))

md.hydrology.steps_per_step = 5
md.smb.steps_per_step = 7  #md.hydrology.steps_per_step
md.timestepping.time_step = 2.
md.timestepping.final_time = 20.0

smb_step = md.timestepping.time_step / md.smb.steps_per_step
duration = np.arange(md.timestepping.start_time, md.timestepping.final_time + smb_step, smb_step)

ddf = 10.0e-3
md.smb.accuref = np.array([[0.5, 0.5], [0., 2.0]])
md.smb.accualti = 0.0
md.smb.accugrad = 0.0

md.smb.runoffref = 0.9 * duration * ddf
md.smb.runoffref = np.vstack((md.smb.runoffref, duration))
md.smb.runoffalti = 0.0
md.smb.runoffgrad = -6.5e-3 * ddf  #lapse rate *ddf*day per year

md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices))
md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices))

md = solve(md, 'Transient')

field_names = ['SedimentWaterHead1', 'EplWaterHead1', 'SedimentHeadResidual1',
               'SedimentWaterHead4', 'EplWaterHead4', 'SedimentHeadResidual4',
               'SedimentWaterHead5', 'EplWaterHead5', 'SedimentHeadResidual5',
               'SedimentWaterHead9', 'EplWaterHead9', 'SedimentHeadResidual9',
               'EplWaterHead10', 'EplWaterHeadSubstep10', 'SedimentWaterHead10',
               'SedimentWaterHeadSubstep10']
field_tolerances = [1e-13, 1e-13, 1e-13,
                    1e-13, 1e-13, 1e-13,
                    1e-13, 5e-12, 1e-11,
                    1e-13, 5e-12, 1e-11,
                    1e-13, 1e-13, 1e-13,
                    1e-13]
field_values = [md.results.TransientSolution[0].SedimentHead,
                md.results.TransientSolution[0].EplHead,
                md.results.TransientSolution[0].SedimentHeadResidual,
                md.results.TransientSolution[3].SedimentHead,
                md.results.TransientSolution[3].EplHead,
                md.results.TransientSolution[3].SedimentHeadResidual,
                md.results.TransientSolution[4].SedimentHead,
                md.results.TransientSolution[4].EplHead,
                md.results.TransientSolution[4].SedimentHeadResidual,
                md.results.TransientSolution[8].SedimentHead,
                md.results.TransientSolution[8].EplHead,
                md.results.TransientSolution[8].SedimentHeadResidual,
                md.results.TransientSolution[-1].EplHead,
                md.results.TransientSolution[-1].EplHeadSubstep,
                md.results.TransientSolution[-1].SedimentHead,
                md.results.TransientSolution[-1].SedimentHeadSubstep]
