#Test Name: ISMIPCFS
import numpy as np
from model import *
from socket import gethostname
from triangle import *
from setmask import *
from parameterize import *
from setflowequation import *
from solve import *

"""
This test is a test from the ISMP-HOM Intercomparison project.
Pattyn and Payne 2006
"""

L_list=[80000.]
results=[]

for L in L_list:
	md=triangle(model(),"../Exp/Square_%d.exp" % L,L/10.)    #size 3*L 
	md=setmask(md,'','')    #ice sheet test
	md=parameterize(md,'../Par/ISMIPC.py')
	md.friction.coefficient=np.sqrt(md.constants.yts*(1000.+1000.*np.sin(md.mesh.x*2.*np.pi/L)*np.sin(md.mesh.y*2.*np.pi/L)))
	md.extrude(10,1.)

	#Add spc on the borders
	pos=np.where(np.logical_or.reduce((md.mesh.x==0.,md.mesh.x==np.max(md.mesh.x),md.mesh.y==0.,md.mesh.y==np.max(md.mesh.y))))
	md.stressbalance.spcvx[pos]=0.
	md.stressbalance.spcvy[pos]=0.
	if   (L==5000.):
		md.stressbalance.spcvx[pos]=15.66
		md.stressbalance.spcvy[pos]=-0.1967
	elif (L==10000.):
		md.stressbalance.spcvx[pos]=16.04
		md.stressbalance.spcvy[pos]=-0.1977
	elif (L==20000.):
		md.stressbalance.spcvx[pos]=16.53
		md.stressbalance.spcvy[pos]=-1.27
	elif (L==40000.):
		md.stressbalance.spcvx[pos]=17.23
		md.stressbalance.spcvy[pos]=-3.17
	elif (L==80000.):
		md.stressbalance.spcvx[pos]=16.68
		md.stressbalance.spcvy[pos]=-2.69
	elif (L==160000.):
		md.stressbalance.spcvx[pos]=16.03
		md.stressbalance.spcvy[pos]=-1.27

	md=setflowequation(md,'FS','all')

	#Compute the stressbalance
	md.cluster=generic('name',gethostname(),'np',8)
	md=solve(md,'Stressbalance')

	#Plot the results and save them
	vx=md.results.StressbalanceSolution.Vx
	vy=md.results.StressbalanceSolution.Vy
	vz=md.results.StressbalanceSolution.Vz
	results.append(md.results.StressbalanceSolution)

#	plotmodel(md,'data',vx,'data',vy,'data',vz,'layer#all',md.mesh.numberoflayers)

#Fields and tolerances to track changes
field_names     =['Vx80km','Vy80km','Vz80km']
field_tolerances=[1e-12,1e-12,1e-12]
field_values=[]
for result in results:
	field_values=field_values+[result.Vx,result.Vy,result.Vz]