Context Navigation

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

Last change on this file since 16560 was 16560, checked in by Mathieu Morlighem, 12 years ago
merged trunk-jpl and trunk for revision 16554
File size: 10.5 KB

Rev	Line
[13020]	1	import numpy
	2	import sys
[13740]	3	import copy
[12038]	4	from fielddisplay import fielddisplay
[13020]	5	from EnumDefinitions import *
	6	from checkfield import *
	7	from WriteData import *
[16560]	8	from MatlabFuncs import *
[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	#}}}
[14640]	77	def setdefaultparameters(self): # {{{
[12123]	78	#maximum of non-linear iterations.
[13020]	79	self.maxiter=100
[12123]	80
	81	#Convergence criterion: absolute, relative and residual
[13020]	82	self.restol=10**-4
	83	self.reltol=0.01
	84	self.abstol=10
[12123]	85
[15564]	86	self.FSreconditioning=10**13
[13020]	87	self.shelf_dampening=0
[12123]	88
	89	#Penalty factor applied kappa=max(stiffness matrix)*10^penalty_factor
[13020]	90	self.penalty_factor=3
[12123]	91
	92	#coefficient to update the viscosity between each iteration of
[15614]	93	#a stressbalance according to the following formula
[12123]	94	#viscosity(n)=viscosity(n)+viscosity_overshoot(viscosity(n)-viscosity(n-1))
[13020]	95	self.viscosity_overshoot=0
[12123]	96
	97	#Stop the iterations of rift if below a threshold
[13020]	98	self.rift_penalty_threshold=0
[12123]	99
	100	#in some solutions, it might be needed to stop a run when only
	101	#a few constraints remain unstable. For thermal computation, this
	102	#parameter is often used.
[13020]	103	self.rift_penalty_lock=10
[12123]	104
[16560]	105	#output default:
	106	self.requested_outputs=['default']
	107
[13020]	108	return self
[12123]	109	#}}}
[16560]	110	def defaultoutputs(self,md): # {{{
	111
	112	if strcmp(md.mesh.meshtype(),'3D'):
	113	list = ['Vx','Vy','Vz','Vel','Pressure']
	114	elif strcmp(md.mesh.meshtype(),'2Dhorizontal'):
	115	list = ['Vx','Vy','Vel','Pressure']
	116	elif strcmp(md.mesh.meshtype(),'2Dvertical'):
	117	list = ['Vx','Vy','Vel','Pressure']
	118	else:
	119	raise TypeError('mesh type not supported yet');
	120	return list
	121
	122	#}}}
[13020]	123	def checkconsistency(self,md,solution,analyses): # {{{
	124
	125	#Early return
[15771]	126	if StressbalanceAnalysisEnum() not in analyses:
[13020]	127	return md
	128
[15771]	129	md = checkfield(md,'stressbalance.spcvx','forcing',1)
	130	md = checkfield(md,'stressbalance.spcvy','forcing',1)
[16560]	131	if strcmp(md.mesh.meshtype(),'3D'):
[15771]	132	md = checkfield(md,'stressbalance.spcvz','forcing',1)
	133	md = checkfield(md,'stressbalance.restol','size',[1],'>',0)
	134	md = checkfield(md,'stressbalance.reltol','size',[1])
	135	md = checkfield(md,'stressbalance.abstol','size',[1])
	136	md = checkfield(md,'stressbalance.isnewton','numel',[1],'values',[0,1,2])
	137	md = checkfield(md,'stressbalance.FSreconditioning','size',[1],'NaN',1)
	138	md = checkfield(md,'stressbalance.viscosity_overshoot','size',[1],'NaN',1)
	139	md = checkfield(md,'stressbalance.maxiter','size',[1],'>=',1)
	140	md = checkfield(md,'stressbalance.referential','size',[md.mesh.numberofvertices,6])
	141	md = checkfield(md,'stressbalance.loadingforce','size',[md.mesh.numberofvertices,3])
[16560]	142	md = checkfield(md,'stressbalance.requested_outputs','stringrow',1);
[13020]	143
	144	#singular solution
[15771]	145	# if ~any((~isnan(md.stressbalance.spcvx)+~isnan(md.stressbalance.spcvy))==2),
	146	if not numpy.any(numpy.logical_and(numpy.logical_not(numpy.isnan(md.stressbalance.spcvx)),numpy.logical_not(numpy.isnan(md.stressbalance.spcvy)))):
[13020]	147	md.checkmessage("model is not well posed (singular). You need at least one node with fixed velocity!")
	148	#CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES
[15771]	149	# if any(sum(isnan(md.stressbalance.referential),2)~=0 & sum(isnan(md.stressbalance.referential),2)~=6),
	150	if numpy.any(numpy.logical_and(numpy.sum(numpy.isnan(md.stressbalance.referential),axis=1)!=0,numpy.sum(numpy.isnan(md.stressbalance.referential),axis=1)!=6)):
	151	md.checkmessage("Each line of stressbalance.referential should contain either only NaN values or no NaN values")
[13020]	152	#CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL
[15771]	153	# if any(sum(isnan(md.stressbalance.referential),2)==0),
	154	if numpy.any(numpy.sum(numpy.isnan(md.stressbalance.referential),axis=1)==0):
	155	pos=[i for i,item in enumerate(numpy.sum(numpy.isnan(md.stressbalance.referential),axis=1)) if item==0]
[13020]	156	# numpy.inner (and numpy.dot) calculate all the dot product permutations, resulting in a full matrix multiply
[15771]	157	# if numpy.any(numpy.abs(numpy.inner(md.stressbalance.referential[pos,0:2],md.stressbalance.referential[pos,3:5]).diagonal())>sys.float_info.epsilon):
	158	# md.checkmessage("Vectors in stressbalance.referential (columns 1 to 3 and 4 to 6) must be orthogonal")
	159	for item in md.stressbalance.referential[pos,:]:
[13020]	160	if numpy.abs(numpy.inner(item[0:2],item[3:5]))>sys.float_info.epsilon:
[15771]	161	md.checkmessage("Vectors in stressbalance.referential (columns 1 to 3 and 4 to 6) must be orthogonal")
[13020]	162	#CHECK THAT NO rotation specified for FS Grounded ice at base
[16560]	163	if strcmp(md.mesh.meshtype(),'3D') and md.flowequation.isFS:
[15987]	164	pos=numpy.nonzero(numpy.logical_and(md.mask.groundedice_levelset,md.mesh.vertexonbed))
[15771]	165	if numpy.any(numpy.logical_not(numpy.isnan(md.stressbalance.referential[pos,:]))):
[13020]	166	md.checkmessage("no referential should be specified for basal vertices of grounded ice")
	167
	168	return md
	169	# }}}
[15131]	170	def marshall(self,md,fid): # {{{
[15621]	171
	172	yts=365.024.03600.0
	173
[15771]	174	WriteData(fid,'object',self,'class','stressbalance','fieldname','spcvx','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1)
	175	WriteData(fid,'object',self,'class','stressbalance','fieldname','spcvy','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1)
	176	WriteData(fid,'object',self,'class','stressbalance','fieldname','spcvz','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1)
	177	WriteData(fid,'object',self,'class','stressbalance','fieldname','restol','format','Double')
	178	WriteData(fid,'object',self,'class','stressbalance','fieldname','reltol','format','Double')
	179	WriteData(fid,'object',self,'class','stressbalance','fieldname','abstol','format','Double')
	180	WriteData(fid,'object',self,'class','stressbalance','fieldname','isnewton','format','Integer')
	181	WriteData(fid,'object',self,'class','stressbalance','fieldname','FSreconditioning','format','Double')
	182	WriteData(fid,'object',self,'class','stressbalance','fieldname','viscosity_overshoot','format','Double')
	183	WriteData(fid,'object',self,'class','stressbalance','fieldname','maxiter','format','Integer')
	184	WriteData(fid,'object',self,'class','stressbalance','fieldname','shelf_dampening','format','Integer')
	185	WriteData(fid,'object',self,'class','stressbalance','fieldname','vertex_pairing','format','DoubleMat','mattype',3)
	186	WriteData(fid,'object',self,'class','stressbalance','fieldname','penalty_factor','format','Double')
	187	WriteData(fid,'object',self,'class','stressbalance','fieldname','rift_penalty_lock','format','Integer')
	188	WriteData(fid,'object',self,'class','stressbalance','fieldname','rift_penalty_threshold','format','Integer')
	189	WriteData(fid,'object',self,'class','stressbalance','fieldname','referential','format','DoubleMat','mattype',1)
[16560]	190
[14529]	191	WriteData(fid,'data',self.loadingforce[:,0],'format','DoubleMat','mattype',1,'enum',LoadingforceXEnum())
	192	WriteData(fid,'data',self.loadingforce[:,1],'format','DoubleMat','mattype',1,'enum',LoadingforceYEnum())
	193	WriteData(fid,'data',self.loadingforce[:,2],'format','DoubleMat','mattype',1,'enum',LoadingforceZEnum())
[16560]	194
	195	#process requested outputs
	196	outputs = self.requested_outputs
	197	indices = [i for i, x in enumerate(outputs) if x == 'default']
	198	if len(indices) > 0:
	199	outputscopy=outputs[0:max(0,indices[0]-1)]+self.defaultoutputs(md)+outputs[indices[0]+1:]
	200	outputs =outputscopy
	201	WriteData(fid,'data',outputs,'enum',StressbalanceRequestedOutputsEnum(),'format','StringArray')
[13020]	202	# }}}

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

Download in other formats: