Context Navigation

source: issm/trunk/src/m/classes/stressbalance.py@ 21341

Last change on this file since 21341 was 21341, checked in by Mathieu Morlighem, 8 years ago
merged trunk-jpl and trunk for revision 21337
File size: 11.3 KB

Rev	Line
[21341]	1	import numpy as np
[13020]	2	import sys
[13740]	3	import copy
[19105]	4	from project3d import project3d
[12038]	5	from fielddisplay import fielddisplay
[17806]	6	from checkfield import checkfield
	7	from WriteData import WriteData
	8	import MatlabFuncs as m
[12038]	9
[15614]	10	class stressbalance(object):
[13020]	11	"""
[15614]	12	STRESSBALANCE class definition
[13020]	13
	14	Usage:
[15614]	15	stressbalance=stressbalance();
[13020]	16	"""
	17
[14640]	18	def __init__(self): # {{{
[12038]	19	self.spcvx = float('NaN')
	20	self.spcvy = float('NaN')
	21	self.spcvz = float('NaN')
	22	self.restol = 0
	23	self.reltol = 0
	24	self.abstol = 0
	25	self.isnewton = 0
[15564]	26	self.FSreconditioning = 0
[12038]	27	self.viscosity_overshoot = 0
	28	self.icefront = float('NaN')
	29	self.maxiter = 0
	30	self.shelf_dampening = 0
	31	self.vertex_pairing = float('NaN')
	32	self.penalty_factor = float('NaN')
	33	self.rift_penalty_lock = float('NaN')
	34	self.rift_penalty_threshold = 0
	35	self.referential = float('NaN')
[14529]	36	self.loadingforce = float('NaN')
[16560]	37	self.requested_outputs = []
[12123]	38
	39	#set defaults
	40	self.setdefaultparameters()
	41
[12038]	42	#}}}
[14640]	43	def __repr__(self): # {{{
[12038]	44
[15614]	45	string=' StressBalance solution parameters:'
[14141]	46	string="%s\n%s"%(string,' Convergence criteria:')
[13020]	47	string="%s\n%s"%(string,fielddisplay(self,'restol','mechanical equilibrium residual convergence criterion'))
[14640]	48	string="%s\n%s"%(string,fielddisplay(self,'reltol','velocity relative convergence criterion, NaN: not applied'))
	49	string="%s\n%s"%(string,fielddisplay(self,'abstol','velocity absolute convergence criterion, NaN: not applied'))
[14413]	50	string="%s\n%s"%(string,fielddisplay(self,'isnewton',"0: Picard's fixed point, 1: Newton's method, 2: hybrid"))
[13020]	51	string="%s\n%s"%(string,fielddisplay(self,'maxiter','maximum number of nonlinear iterations'))
	52	string="%s\n%s"%(string,fielddisplay(self,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)'))
[12038]	53
[14141]	54	string="%s\n%s"%(string,'\n boundary conditions:')
[14640]	55	string="%s\n%s"%(string,fielddisplay(self,'spcvx','x-axis velocity constraint (NaN means no constraint) [m/yr]'))
	56	string="%s\n%s"%(string,fielddisplay(self,'spcvy','y-axis velocity constraint (NaN means no constraint) [m/yr]'))
	57	string="%s\n%s"%(string,fielddisplay(self,'spcvz','z-axis velocity constraint (NaN means no constraint) [m/yr]'))
	58	string="%s\n%s"%(string,fielddisplay(self,'icefront','segments on ice front list (last column 0: Air, 1: Water, 2: Ice'))
[12038]	59
[14141]	60	string="%s\n%s"%(string,'\n Rift options:')
[13020]	61	string="%s\n%s"%(string,fielddisplay(self,'rift_penalty_threshold','threshold for instability of mechanical constraints'))
	62	string="%s\n%s"%(string,fielddisplay(self,'rift_penalty_lock','number of iterations before rift penalties are locked'))
[12038]	63
[14141]	64	string="%s\n%s"%(string,'\n Penalty options:')
[13020]	65	string="%s\n%s"%(string,fielddisplay(self,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset'))
	66	string="%s\n%s"%(string,fielddisplay(self,'vertex_pairing','pairs of vertices that are penalized'))
[12038]	67
[14141]	68	string="%s\n%s"%(string,'\n Other:')
[15564]	69	string="%s\n%s"%(string,fielddisplay(self,'shelf_dampening','use dampening for floating ice ? Only for FS model'))
	70	string="%s\n%s"%(string,fielddisplay(self,'FSreconditioning','multiplier for incompressibility equation. Only for FS model'))
[13020]	71	string="%s\n%s"%(string,fielddisplay(self,'referential','local referential'))
[14640]	72	string="%s\n%s"%(string,fielddisplay(self,'loadingforce','loading force applied on each point [N/m^3]'))
[13020]	73	string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested'))
[12038]	74
	75	return string
	76	#}}}
[19105]	77	def extrude(self,md): # {{{
	78	self.spcvx=project3d(md,'vector',self.spcvx,'type','node')
	79	self.spcvy=project3d(md,'vector',self.spcvy,'type','node')
	80	self.spcvz=project3d(md,'vector',self.spcvz,'type','node')
	81	self.referential=project3d(md,'vector',self.referential,'type','node')
	82	self.loadingforce=project3d(md,'vector',self.loadingforce,'type','node')
	83
	84	return self
	85	#}}}
[14640]	86	def setdefaultparameters(self): # {{{
[12123]	87	#maximum of non-linear iterations.
[13020]	88	self.maxiter=100
[12123]	89
	90	#Convergence criterion: absolute, relative and residual
[13020]	91	self.restol=10**-4
	92	self.reltol=0.01
	93	self.abstol=10
[12123]	94
[15564]	95	self.FSreconditioning=10**13
[13020]	96	self.shelf_dampening=0
[12123]	97
	98	#Penalty factor applied kappa=max(stiffness matrix)*10^penalty_factor
[13020]	99	self.penalty_factor=3
[12123]	100
	101	#coefficient to update the viscosity between each iteration of
[15614]	102	#a stressbalance according to the following formula
[12123]	103	#viscosity(n)=viscosity(n)+viscosity_overshoot(viscosity(n)-viscosity(n-1))
[13020]	104	self.viscosity_overshoot=0
[12123]	105
	106	#Stop the iterations of rift if below a threshold
[13020]	107	self.rift_penalty_threshold=0
[12123]	108
	109	#in some solutions, it might be needed to stop a run when only
	110	#a few constraints remain unstable. For thermal computation, this
	111	#parameter is often used.
[13020]	112	self.rift_penalty_lock=10
[12123]	113
[16560]	114	#output default:
	115	self.requested_outputs=['default']
	116
[13020]	117	return self
[12123]	118	#}}}
[16560]	119	def defaultoutputs(self,md): # {{{
	120
[17806]	121	if md.mesh.dimension()==3:
[16560]	122	list = ['Vx','Vy','Vz','Vel','Pressure']
[17806]	123	else:
[16560]	124	list = ['Vx','Vy','Vel','Pressure']
	125	return list
	126
	127	#}}}
[13020]	128	def checkconsistency(self,md,solution,analyses): # {{{
	129
	130	#Early return
[21341]	131	if 'StressbalanceAnalysis' not in analyses:
[13020]	132	return md
	133
[20500]	134	md = checkfield(md,'fieldname','stressbalance.spcvx','Inf',1,'timeseries',1)
	135	md = checkfield(md,'fieldname','stressbalance.spcvy','Inf',1,'timeseries',1)
[17806]	136	if m.strcmp(md.mesh.domaintype(),'3D'):
[20500]	137	md = checkfield(md,'fieldname','stressbalance.spcvz','Inf',1,'timeseries',1)
[17806]	138	md = checkfield(md,'fieldname','stressbalance.restol','size',[1],'>',0)
	139	md = checkfield(md,'fieldname','stressbalance.reltol','size',[1])
	140	md = checkfield(md,'fieldname','stressbalance.abstol','size',[1])
	141	md = checkfield(md,'fieldname','stressbalance.isnewton','numel',[1],'values',[0,1,2])
[20500]	142	md = checkfield(md,'fieldname','stressbalance.FSreconditioning','size',[1],'NaN',1,'Inf',1)
	143	md = checkfield(md,'fieldname','stressbalance.viscosity_overshoot','size',[1],'NaN',1,'Inf',1)
[17806]	144	md = checkfield(md,'fieldname','stressbalance.maxiter','size',[1],'>=',1)
	145	md = checkfield(md,'fieldname','stressbalance.referential','size',[md.mesh.numberofvertices,6])
	146	md = checkfield(md,'fieldname','stressbalance.loadingforce','size',[md.mesh.numberofvertices,3])
	147	md = checkfield(md,'fieldname','stressbalance.requested_outputs','stringrow',1);
[13020]	148
	149	#singular solution
[15771]	150	# if ~any((~isnan(md.stressbalance.spcvx)+~isnan(md.stressbalance.spcvy))==2),
[21341]	151	if not np.any(np.logical_and(np.logical_not(np.isnan(md.stressbalance.spcvx)),np.logical_not(np.isnan(md.stressbalance.spcvy)))):
[20500]	152	print "\n !!! Warning: no spc applied, model might not be well posed if no basal friction is applied, check for solution crash\n"
[13020]	153	#CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES
[15771]	154	# if any(sum(isnan(md.stressbalance.referential),2)~=0 & sum(isnan(md.stressbalance.referential),2)~=6),
[21341]	155	if np.any(np.logical_and(np.sum(np.isnan(md.stressbalance.referential),axis=1)!=0,np.sum(np.isnan(md.stressbalance.referential),axis=1)!=6)):
[15771]	156	md.checkmessage("Each line of stressbalance.referential should contain either only NaN values or no NaN values")
[13020]	157	#CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL
[15771]	158	# if any(sum(isnan(md.stressbalance.referential),2)==0),
[21341]	159	if np.any(np.sum(np.isnan(md.stressbalance.referential),axis=1)==0):
	160	pos=[i for i,item in enumerate(np.sum(np.isnan(md.stressbalance.referential),axis=1)) if item==0]
	161	# np.inner (and np.dot) calculate all the dot product permutations, resulting in a full matrix multiply
	162	# if np.any(np.abs(np.inner(md.stressbalance.referential[pos,0:2],md.stressbalance.referential[pos,3:5]).diagonal())>sys.float_info.epsilon):
[15771]	163	# md.checkmessage("Vectors in stressbalance.referential (columns 1 to 3 and 4 to 6) must be orthogonal")
	164	for item in md.stressbalance.referential[pos,:]:
[21341]	165	if np.abs(np.inner(item[0:2],item[3:5]))>sys.float_info.epsilon:
[15771]	166	md.checkmessage("Vectors in stressbalance.referential (columns 1 to 3 and 4 to 6) must be orthogonal")
[13020]	167	#CHECK THAT NO rotation specified for FS Grounded ice at base
[17806]	168	if m.strcmp(md.mesh.domaintype(),'3D') and md.flowequation.isFS:
[21341]	169	pos=np.nonzero(np.logical_and(md.mask.groundedice_levelset,md.mesh.vertexonbase))
	170	if np.any(np.logical_not(np.isnan(md.stressbalance.referential[pos,:]))):
[13020]	171	md.checkmessage("no referential should be specified for basal vertices of grounded ice")
	172
	173	return md
	174	# }}}
[21341]	175	def marshall(self,prefix,md,fid): # {{{
[15621]	176
[21341]	177	WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','vertex_pairing','format','DoubleMat','mattype',3)
[15621]	178
[21341]	179	yts=md.constants.yts
	180
	181	WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','spcvx','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
	182	WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','spcvy','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
	183	WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','spcvz','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
	184	WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','restol','format','Double')
	185	WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','reltol','format','Double')
	186	WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','abstol','format','Double','scale',1./yts)
	187	WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','isnewton','format','Integer')
	188	WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','FSreconditioning','format','Double')
	189	WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','viscosity_overshoot','format','Double')
	190	WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','maxiter','format','Integer')
	191	WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','shelf_dampening','format','Integer')
	192	WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','penalty_factor','format','Double')
	193	WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','rift_penalty_lock','format','Integer')
	194	WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','rift_penalty_threshold','format','Integer')
	195	WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','referential','format','DoubleMat','mattype',1)
[20500]	196
[21341]	197	if isinstance(self.loadingforce, (list, tuple, np.ndarray)):
[20500]	198	lx=self.loadingforce[:,0];
	199	ly=self.loadingforce[:,1];
	200	lz=self.loadingforce[:,2];
	201	else:
	202	lx=float('NaN'); ly=float('NaN'); lz=float('NaN');
[16560]	203
[21341]	204	WriteData(fid,prefix,'data',lx,'format','DoubleMat','mattype',1,'name','md.stressbalance.loadingforcex')
	205	WriteData(fid,prefix,'data',ly,'format','DoubleMat','mattype',1,'name','md.stressbalance.loadingforcey')
	206	WriteData(fid,prefix,'data',lz,'format','DoubleMat','mattype',1,'name','md.stressbalance.loadingforcez')
[16560]	207
	208	#process requested outputs
	209	outputs = self.requested_outputs
	210	indices = [i for i, x in enumerate(outputs) if x == 'default']
	211	if len(indices) > 0:
	212	outputscopy=outputs[0:max(0,indices[0]-1)]+self.defaultoutputs(md)+outputs[indices[0]+1:]
	213	outputs =outputscopy
[21341]	214	WriteData(fid,prefix,'data',outputs,'name','md.stressbalance.requested_outputs','format','StringArray')
[13020]	215	# }}}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: