Context Navigation

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

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

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

Download in other formats: