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

Last change on this file since 26744 was 26744, checked in by Mathieu Morlighem, 3 years ago

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