source: issm/oecreview/Archive/24684-25833/ISSM-25302-25303.diff

Last change on this file was 25834, checked in by Mathieu Morlighem, 4 years ago

CHG: added 24684-25833
File size: 20.7 KB
RevLine 
[25834]1Index: ../trunk-jpl/test/NightlyRun/test2085.m
2===================================================================
3--- ../trunk-jpl/test/NightlyRun/test2085.m (revision 25302)
4+++ ../trunk-jpl/test/NightlyRun/test2085.m (revision 25303)
5@@ -10,10 +10,10 @@
6 md.materials.numlayers = 10;
7 md.love.forcing_type = 9;
8
9- md.materials.density= (1:md.materials.numlayers)'*0+5511;
10- md.materials.lame_mu= (1:md.materials.numlayers)'*0+0.75e11;
11- md.materials.viscosity=(1:md.materials.numlayers)'*0+1e21;
12- md.materials.lame_lambda=md.materials.lame_mu*0+5e17;
13+ md.materials.density=zeros(md.materials.numlayers,1)+5511;
14+ md.materials.lame_mu=zeros(md.materials.numlayers,1)+0.75e11;
15+ md.materials.viscosity=zeros(md.materials.numlayers,1)+1e21;
16+ md.materials.lame_lambda=zeros(md.materials.numlayers,1)+5e17;
17 md.materials.issolid=ones(md.materials.numlayers,1);
18 md.materials.isburgers=zeros(md.materials.numlayers,1);
19 md.materials.burgers_mu=md.materials.lame_mu/3;
20@@ -24,7 +24,7 @@
21
22 md.love.allow_layer_deletion=1;
23 md.love.frequencies=0;
24- md.love.nfreq=length(md.love.frequencies);
25+ md.love.nfreq=length(md.love.frequencies); % TODO: Why are we grabbing length of a scalar here?
26
27 md.love.sh_nmin = 2;
28 md.love.sh_nmax = 2;
29@@ -37,7 +37,7 @@
30 md=solve(md,'lv');
31
32 % extract love kernels {{{
33- y1=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,1)));
34+ y1=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,1)));
35 y2=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,2)));
36 y3=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,3)));
37 y4=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,4)));
38@@ -147,18 +147,18 @@
39 'y6_load_degree2_interface1','y6_load_degree2_interface2','y6_load_degree2_interface7','y6_load_degree2_interface10','y6_load_degree2_interface11',...
40 };
41 field_tolerances={...
42- 1e-10,1e-10,1e-10,1e-10,1e-10,...
43- 1e-10,1e-10,1e-10,1e-10,1e-10,...
44- 1e-10,1e-10,1e-10,1e-10,1e-10,...
45- 1e-10,1e-10,1e-10,1e-10,1e-10,...
46- 1e-10,1e-10,1e-10,1e-10,1e-10,...
47- 1e-10,1e-10,1e-10,1e-10,1e-10,...
48- 1e-10,1e-10,1e-10,1e-10,1e-10,...
49- 1e-10,1e-10,1e-10,1e-10,1e-10,...
50- 1e-10,1e-10,1e-10,1e-10,1e-10,...
51- 1e-10,1e-10,1e-10,1e-10,1e-10,...
52- 1e-10,1e-10,1e-10,1e-10,1e-10,...
53- 1e-10,1e-10,1e-10,1e-10,1e-10,...
54+ 9e-8, 1e-7, 3e-7, 3e-7, 3e-7,...
55+ 9e-8, 1e-7, 3e-7, 3e-7, 1e-10,...
56+ 9e-8, 8e-8, 2e-8, 2e-7, 4e-7,...
57+ 9e-8, 9e-8, 2e-7, 4e-7, 1e-10,...
58+ 4e-7, 4e-7, 2e-7, 3e-8, 2e-8,...
59+ 2e-5, 2e-6, 2e-6, 1e-6, 2e-7,...
60+ 3e-6, 3e-6, 3e-6, 4e-6, 4e-6,...
61+ 3e-6, 3e-6, 2e-6, 6e-7, 1e-10,...
62+ 3e-6, 3e-6, 5e-7, 3e-6, 5e-6,...
63+ 3e-6, 3e-6, 9e-7, 7e-7, 1e-10,...
64+ 4e-6, 4e-6, 3e-6, 5e-7, 3e-7,...
65+ 3e-6, 3e-6, 2e-6, 7e-7, 2e-7...
66 };
67 field_values={...
68 y1_tidal_degree2_interface1, y1_tidal_degree2_interface2, y1_tidal_degree2_interface7, y1_tidal_degree2_interface10,y1_tidal_degree2_interface11,...
69Index: ../trunk-jpl/test/NightlyRun/test2085.py
70===================================================================
71--- ../trunk-jpl/test/NightlyRun/test2085.py (revision 25302)
72+++ ../trunk-jpl/test/NightlyRun/test2085.py (revision 25303)
73@@ -1,49 +1,188 @@
74-#Test Name: LovenumberstAtDepth.
75-#Same as test #1 of test2084.m
76+# Test Name: FourierLoveKernels
77+# Homogenous Earth, for which analytic solutions exist.
78+# Love kernels for degree 2 (tested against analytic solns).
79
80-from model import *
81+
82 from socket import gethostname
83-from solve import *
84-from numpy import *
85+
86+import numpy as np
87+
88 from generic import generic
89 from materials import *
90+from model import *
91+from solve import *
92
93+# Fro volumetric potential
94 md = model()
95-md.cluster = generic('name', gethostname(), 'np', 1)
96
97 md.materials = materials('litho')
98-md.miscellaneous.name = 'FourierLoveTest'
99-md.groundingline.migration = 'None'
100
101-md.verbose = verbose('111111101')
102-cst = 365.25 * 24 * 3600 * 1000
103+md.materials.numlayers = 10
104+md.love.forcing_type = 9
105
106-md.materials.numlayers = 6
107-md.materials.radius = np.array([10, 1222.5, 3.4800e+03, 5.7010e+03, 5.9510e+03, 6.3010e+03, 6.3710e+03]).reshape(-1, 1) * 1e3
108-md.materials.density = np.array([1.0750e4, 1.0750e+04, 4.9780e+03, 3.8710e+03, 3.4380e+03, 3.0370e+03]).reshape(-1, 1)
109-md.materials.lame_mu = np.array([1e-5, 0, 2.2834e+00, 1.0549e+00, 7.0363e-01, 5.0605e-01]).reshape(-1, 1) * 1e11
110-md.materials.viscosity = np.array([0, 0, 2.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+25]).reshape(-1, 1) * 1e21
111-md.materials.lame_lambda = np.array(md.materials.lame_mu) * 0 + 5e14
112-md.materials.issolid = np.array([1, 0, 1, 1, 1, 1]).reshape(-1, 1)
113-md.materials.isburgers = np.zeros((md.materials.numlayers)).reshape(-1, 1)
114+md.materials.density = np.zeros((md.materials.numlayers, 1)) + 5511
115+md.materials.lame_mu = np.zeros((md.materials.numlayers, 1)) + 0.75e11
116+md.materials.viscosity = np.zeros((md.materials.numlayers, 1)) + 1e21
117+md.materials.lame_lambda = np.zeros((md.materials.numlayers, 1)) + 0.5e17
118+md.materials.issolid = np.ones((md.materials.numlayers, 1))
119+md.materials.isburgers = np.zeros((md.materials.numlayers, 1))
120+md.materials.burgers_mu = md.materials.lame_mu / 3
121+md.materials.burgers_viscosity = md.materials.viscosity / 10
122
123-md.love.love_kernels = 1
124+md.materials.radius = np.linspace(10e3, 6371e3, md.materials.numlayers + 1).reshape(-1, 1)
125+md.love.g0 = 9.8134357285509388 # directly grabbed from fourierlovesolver for this particular case
126+
127 md.love.allow_layer_deletion = 1
128-md.love.frequencies = (np.array([0]) * 2 * pi).reshape(-1, 1) / cst
129-md.love.nfreq = len(md.love.frequencies)
130+md.love.frequencies = 0
131+md.love.nfreq = 1
132+
133+md.love.sh_nmin = 2
134 md.love.sh_nmax = 2
135+md.love.love_kernels = 1
136
137-md.materials.burgers_mu = md.materials.lame_mu
138-md.materials.burgers_viscosity = md.materials.viscosity
139+md.miscellaneous.name = 'kernels'
140+md.cluster = generic('name', gethostname(), 'np', 1)
141+md.verbose = verbose('111111101')
142
143-md = solve(md, 'lv')
144+md = solve(md,'lv')
145
146+# Extract love kernels #{{{
147+y1 = md.results.LoveSolution.LoveKernelsReal[:,0,:,0].squeeze()
148+y2 = md.results.LoveSolution.LoveKernelsReal[:,0,:,1].squeeze()
149+y3 = md.results.LoveSolution.LoveKernelsReal[:,0,:,2].squeeze()
150+y4 = md.results.LoveSolution.LoveKernelsReal[:,0,:,3].squeeze()
151+y5 = md.results.LoveSolution.LoveKernelsReal[:,0,:,4].squeeze()
152+y6 = md.results.LoveSolution.LoveKernelsReal[:,0,:,5].squeeze()
153+
154+y1_tidal_degree2_interface1 = y1[2, 0]
155+y1_tidal_degree2_interface2 = y1[2, 1]
156+y1_tidal_degree2_interface7 = y1[2, 6]
157+y1_tidal_degree2_interface10 = y1[2, 9]
158+y1_tidal_degree2_interface11 = y1[2, 10]
159+
160+y2_tidal_degree2_interface1 = y2[2, 0]
161+y2_tidal_degree2_interface2 = y2[2, 1]
162+y2_tidal_degree2_interface7 = y2[2, 6]
163+y2_tidal_degree2_interface10 = y2[2, 9]
164+y2_tidal_degree2_interface11 = y2[2, 10]
165+
166+y3_tidal_degree2_interface1 = y3[2, 0]
167+y3_tidal_degree2_interface2 = y3[2, 1]
168+y3_tidal_degree2_interface7 = y3[2, 6]
169+y3_tidal_degree2_interface10 = y3[2, 9]
170+y3_tidal_degree2_interface11 = y3[2, 10]
171+
172+y4_tidal_degree2_interface1 = y4[2, 0]
173+y4_tidal_degree2_interface2 = y4[2, 1]
174+y4_tidal_degree2_interface7 = y4[2, 6]
175+y4_tidal_degree2_interface10 = y4[2, 9]
176+y4_tidal_degree2_interface11 = y4[2, 10]
177+
178+y5_tidal_degree2_interface1 = y5[2, 0]
179+y5_tidal_degree2_interface2 = y5[2, 1]
180+y5_tidal_degree2_interface7 = y5[2, 6]
181+y5_tidal_degree2_interface10 = y5[2, 9]
182+y5_tidal_degree2_interface11 = y5[2, 10]
183+
184+y6_tidal_degree2_interface1 = y6[2, 0]
185+y6_tidal_degree2_interface2 = y6[2, 1]
186+y6_tidal_degree2_interface7 = y6[2, 6]
187+y6_tidal_degree2_interface10 = y6[2, 9]
188+y6_tidal_degree2_interface11 = y6[2, 10]
189+#}}}
190+
191+# For surface load
192+md.love.forcing_type = 11
193+
194+md = solve(md,'lv')
195+
196+# Extract love kernels #{{{
197+y1 = md.results.LoveSolution.LoveKernelsReal[:,0,:,0].squeeze()
198+y2 = md.results.LoveSolution.LoveKernelsReal[:,0,:,1].squeeze()
199+y3 = md.results.LoveSolution.LoveKernelsReal[:,0,:,2].squeeze()
200+y4 = md.results.LoveSolution.LoveKernelsReal[:,0,:,3].squeeze()
201+y5 = md.results.LoveSolution.LoveKernelsReal[:,0,:,4].squeeze()
202+y6 = md.results.LoveSolution.LoveKernelsReal[:,0,:,5].squeeze()
203+
204+y1_load_degree2_interface1 = y1[2, 0]
205+y1_load_degree2_interface2 = y1[2, 1]
206+y1_load_degree2_interface7 = y1[2, 6]
207+y1_load_degree2_interface10 = y1[2, 9]
208+y1_load_degree2_interface11 = y1[2, 10]
209+
210+y2_load_degree2_interface1 = y2[2, 0]
211+y2_load_degree2_interface2 = y2[2, 1]
212+y2_load_degree2_interface7 = y2[2, 6]
213+y2_load_degree2_interface10 = y2[2, 9]
214+y2_load_degree2_interface11 = y2[2, 10]
215+
216+y3_load_degree2_interface1 = y3[2, 0]
217+y3_load_degree2_interface2 = y3[2, 1]
218+y3_load_degree2_interface7 = y3[2, 6]
219+y3_load_degree2_interface10 = y3[2, 9]
220+y3_load_degree2_interface11 = y3[2, 10]
221+
222+y4_load_degree2_interface1 = y4[2, 0]
223+y4_load_degree2_interface2 = y4[2, 1]
224+y4_load_degree2_interface7 = y4[2, 6]
225+y4_load_degree2_interface10 = y4[2, 9]
226+y4_load_degree2_interface11 = y4[2, 10]
227+
228+y5_load_degree2_interface1 = y5[2, 0]
229+y5_load_degree2_interface2 = y5[2, 1]
230+y5_load_degree2_interface7 = y5[2, 6]
231+y5_load_degree2_interface10 = y5[2, 9]
232+y5_load_degree2_interface11 = y5[2, 10]
233+
234+y6_load_degree2_interface1 = y6[2, 0]
235+y6_load_degree2_interface2 = y6[2, 1]
236+y6_load_degree2_interface7 = y6[2, 6]
237+y6_load_degree2_interface10 = y6[2, 9]
238+y6_load_degree2_interface11 = y6[2, 10]
239+#}}}
240+
241+
242 #Fields and tolerances to track changes
243 #loading love numbers
244-field_names = ['LoveH_loading_elastic', 'LoveK_loading_elastic', 'LoveL_loading_elastic', 'LoveKernels_degree1', 'LoveKernels_degree2']
245-field_tolerances = [1e-10, 1e-10, 1e-10, 1e-10, 1e-10]
246-field_values = [md.results.LoveSolution.LoveHr[:, 0],
247- md.results.LoveSolution.LoveKr[:, 0],
248- md.results.LoveSolution.LoveLr[:, 0],
249- md.results.LoveSolution.LoveKernelsReal[1, 0, :, :],
250- md.results.LoveSolution.LoveKernelsReal[2, 0, :, :]]
251+field_names = [
252+ 'y1_tidal_degree2_interface1', 'y1_tidal_degree2_interface2', 'y1_tidal_degree2_interface7', 'y1_tidal_degree2_interface10', 'y1_tidal_degree2_interface11',
253+ 'y2_tidal_degree2_interface1', 'y2_tidal_degree2_interface2', 'y2_tidal_degree2_interface7', 'y2_tidal_degree2_interface10', 'y2_tidal_degree2_interface11',
254+ 'y3_tidal_degree2_interface1', 'y3_tidal_degree2_interface2', 'y3_tidal_degree2_interface7', 'y3_tidal_degree2_interface10', 'y3_tidal_degree2_interface11',
255+ 'y4_tidal_degree2_interface1', 'y4_tidal_degree2_interface2', 'y4_tidal_degree2_interface7', 'y4_tidal_degree2_interface10', 'y4_tidal_degree2_interface11',
256+ 'y5_tidal_degree2_interface1', 'y5_tidal_degree2_interface2', 'y5_tidal_degree2_interface7', 'y5_tidal_degree2_interface10', 'y5_tidal_degree2_interface11',
257+ 'y6_tidal_degree2_interface1', 'y6_tidal_degree2_interface2', 'y6_tidal_degree2_interface7', 'y6_tidal_degree2_interface10', 'y6_tidal_degree2_interface11',
258+ 'y1_load_degree2_interface1', 'y1_load_degree2_interface2', 'y1_load_degree2_interface7', 'y1_load_degree2_interface10', 'y1_load_degree2_interface11',
259+ 'y2_load_degree2_interface1', 'y2_load_degree2_interface2', 'y2_load_degree2_interface7', 'y2_load_degree2_interface10', 'y2_load_degree2_interface11',
260+ 'y3_load_degree2_interface1', 'y3_load_degree2_interface2', 'y3_load_degree2_interface7', 'y3_load_degree2_interface10', 'y3_load_degree2_interface11',
261+ 'y4_load_degree2_interface1', 'y4_load_degree2_interface2', 'y4_load_degree2_interface7', 'y4_load_degree2_interface10', 'y4_load_degree2_interface11',
262+ 'y5_load_degree2_interface1', 'y5_load_degree2_interface2', 'y5_load_degree2_interface7', 'y5_load_degree2_interface10', 'y5_load_degree2_interface11',
263+ 'y6_load_degree2_interface1', 'y6_load_degree2_interface2', 'y6_load_degree2_interface7', 'y6_load_degree2_interface10', 'y6_load_degree2_interface11'
264+ ]
265+field_tolerances = [
266+ 9e-8, 1e-7, 3e-7, 3e-7, 3e-7,
267+ 9e-8, 1e-7, 3e-7, 3e-7, 1e-10,
268+ 9e-8, 8e-8, 2e-8, 2e-7, 4e-7,
269+ 9e-8, 9e-8, 2e-7, 4e-7, 1e-10,
270+ 4e-7, 4e-7, 2e-7, 3e-8, 2e-8,
271+ 2e-5, 2e-6, 2e-6, 1e-6, 2e-7,
272+ 3e-6, 3e-6, 3e-6, 4e-6, 4e-6,
273+ 3e-6, 3e-6, 2e-6, 6e-7, 1e-10,
274+ 3e-6, 3e-6, 5e-7, 3e-6, 5e-6,
275+ 3e-6, 3e-6, 9e-7, 7e-7, 1e-10,
276+ 4e-6, 4e-6, 3e-6, 5e-7, 3e-7,
277+ 3e-6, 3e-6, 2e-6, 7e-7, 2e-7
278+ ]
279+field_values = [
280+ y1_tidal_degree2_interface1, y1_tidal_degree2_interface2, y1_tidal_degree2_interface7, y1_tidal_degree2_interface10,y1_tidal_degree2_interface11,
281+ y2_tidal_degree2_interface1, y2_tidal_degree2_interface2, y2_tidal_degree2_interface7, y2_tidal_degree2_interface10,y2_tidal_degree2_interface11,
282+ y3_tidal_degree2_interface1, y3_tidal_degree2_interface2, y3_tidal_degree2_interface7, y3_tidal_degree2_interface10,y3_tidal_degree2_interface11,
283+ y4_tidal_degree2_interface1, y4_tidal_degree2_interface2, y4_tidal_degree2_interface7, y4_tidal_degree2_interface10,y4_tidal_degree2_interface11,
284+ y5_tidal_degree2_interface1, y5_tidal_degree2_interface2, y5_tidal_degree2_interface7, y5_tidal_degree2_interface10,y5_tidal_degree2_interface11,
285+ y6_tidal_degree2_interface1, y6_tidal_degree2_interface2, y6_tidal_degree2_interface7, y6_tidal_degree2_interface10,y6_tidal_degree2_interface11,
286+ y1_load_degree2_interface1, y1_load_degree2_interface2, y1_load_degree2_interface7, y1_load_degree2_interface10,y1_load_degree2_interface11,
287+ y2_load_degree2_interface1, y2_load_degree2_interface2, y2_load_degree2_interface7, y2_load_degree2_interface10,y2_load_degree2_interface11,
288+ y3_load_degree2_interface1, y3_load_degree2_interface2, y3_load_degree2_interface7, y3_load_degree2_interface10,y3_load_degree2_interface11,
289+ y4_load_degree2_interface1, y4_load_degree2_interface2, y4_load_degree2_interface7, y4_load_degree2_interface10,y4_load_degree2_interface11,
290+ y5_load_degree2_interface1, y5_load_degree2_interface2, y5_load_degree2_interface7, y5_load_degree2_interface10,y5_load_degree2_interface11,
291+ y6_load_degree2_interface1, y6_load_degree2_interface2, y6_load_degree2_interface7, y6_load_degree2_interface10,y6_load_degree2_interface11
292+ ]
293Index: ../trunk-jpl/test/NightlyRun/test2052.m
294===================================================================
295--- ../trunk-jpl/test/NightlyRun/test2052.m (revision 25302)
296+++ ../trunk-jpl/test/NightlyRun/test2052.m (revision 25303)
297@@ -30,8 +30,6 @@
298 %% solve for GIA deflection
299 md=solve(md,'Gia');
300
301-pause
302-
303 %Test Name: GiaIvinsBenchmarksAB2dC1
304 U_AB2dC1 = md.results.GiaSolution.UGia;
305 URate_AB2dC1 = md.results.GiaSolution.UGiaRate;
306Index: ../trunk-jpl/src/m/classes/geometry.m
307===================================================================
308--- ../trunk-jpl/src/m/classes/geometry.m (revision 25302)
309+++ ../trunk-jpl/src/m/classes/geometry.m (revision 25303)
310@@ -89,11 +89,8 @@
311
312 end % }}}
313 function marshall(self,prefix,md,fid) % {{{
314- disp(md.geometry.thickness)
315 WriteData(fid,prefix,'object',self,'fieldname','surface','format','DoubleMat','mattype',1);
316 WriteData(fid,prefix,'object',self,'fieldname','thickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
317- disp(md.geometry.thickness)
318- pause
319 WriteData(fid,prefix,'object',self,'fieldname','base','format','DoubleMat','mattype',1);
320 WriteData(fid,prefix,'object',self,'fieldname','bed','format','DoubleMat','mattype',1);
321 WriteData(fid,prefix,'object',self,'fieldname','hydrostatic_ratio','format','DoubleMat','mattype',1);
322Index: ../trunk-jpl/src/m/solve/parseresultsfromdisk.py
323===================================================================
324--- ../trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 25302)
325+++ ../trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 25303)
326@@ -247,34 +247,34 @@
327 mu0 = md.love.mu0
328 rr = md.materials.radius
329 rho = md.materials.density
330- rho_avg = (rho * np.diff(np.power(rr, 3), n=1, axis=0)) / np.diff(np.power(rr, 3).sum()).sum()
331+ rho_avg_partial = np.diff(np.power(rr, 3), n=1, axis=0)
332+ rho_avg = ((rho * rho_avg_partial) / rho_avg_partial.sum()).sum()
333 temp_field = np.empty((degmax + 1, nfreq, nlayer + 1, 6))
334 temp_field.fill(0.0)
335 for ii in range(degmax + 1):
336 for jj in range(nfreq):
337 for kk in range(nlayer + 1):
338- if kk < nlayer + 1:
339+ if kk < nlayer: # NOTE: Upper bound of range is non-inclusive (compare to src/m/solve/parseresultsfromdisk.m)
340 ll = ii * (nlayer + 1) * 6 + (kk * 6 + 1) + 3
341- temp_field[ii, jj, kk, 0] = field[ll + (1 - 1), jj] * r0 # mm = 4
342- temp_field[ii, jj, kk, 1] = field[ll + (2 - 1), jj] * mu0 # mm = 5
343- temp_field[ii, jj, kk, 2] = field[ll + (3 - 1), jj] * r0 # mm = 6
344- temp_field[ii, jj, kk, 3] = field[ll + (4 - 1), jj] * mu0 # mm = 1
345- temp_field[ii, jj, kk, 4] = field[ll + (5 - 1), jj] * r0 * g0 # mm = 2
346- temp_field[ii, jj, kk, 5] = field[ll + (6 - 1), jj] * g0 # mm = 3
347- print(temp_field)
348+ temp_field[ii, jj, kk, 0] = field[ll + (0 - 1), jj] * r0 # mm = 4
349+ temp_field[ii, jj, kk, 1] = field[ll + (1 - 1), jj] * mu0 # mm = 5
350+ temp_field[ii, jj, kk, 2] = field[ll + (2 - 1), jj] * r0 # mm = 6
351+ temp_field[ii, jj, kk, 3] = field[ll + (3 - 1), jj] * mu0 # mm = 1
352+ temp_field[ii, jj, kk, 4] = field[ll + (4 - 1), jj] * r0 * g0 # mm = 2
353+ temp_field[ii, jj, kk, 5] = field[ll + (5 - 1), jj] * g0 # mm = 3
354 else: # surface
355- ll = ii * (nlayer + 1) * 6 - 2
356- temp_field[ii, jj, kk, 0] = field[ll + (1 - 1), jj] * r0
357- temp_field[ii, jj, kk, 2] = field[ll + (2 - 1), jj] * r0
358- temp_field[ii, jj, kk, 4] = field[ll + (3 - 1), jj] * r0 * g0
359+ ll = (ii + 1) * (nlayer + 1) * 6 - 2
360+ temp_field[ii, jj, kk, 0] = field[ll + (0 - 1), jj] * r0
361+ temp_field[ii, jj, kk, 2] = field[ll + (1 - 1), jj] * r0
362+ temp_field[ii, jj, kk, 4] = field[ll + (2 - 1), jj] * r0 * g0
363 # surface BC
364 temp_field[ii, jj, kk, 3] = 0
365 if md.love.forcing_type == 9:
366 temp_field[ii, jj, kk, 1] = 0
367- temp_field[ii, jj, kk, 5] = (2 * ii - 1) / r0 - ii * field[ll + (3 - 1), jj] * g0
368+ temp_field[ii, jj, kk, 5] = (2 * (ii + 1) - 1) / r0 - (ii + 1) * field[ll + (2 - 1), jj] * g0
369 elif md.love.forcing_type == 11:
370- temp_field[ii, jj, kk, 1] = -(2 * (ii - 1) + 1) * rho_avg / 3
371- temp_field[ii, jj, kk, 5] = (2 * ii - 1) / r0 - ii * field[ll + (3 - 1), jj] * g0
372+ temp_field[ii, jj, kk, 1] = -(2 * ii + 1) * rho_avg / 3
373+ temp_field[ii, jj, kk, 5] = (2 * (ii + 1) - 1) / r0 - (ii + 1) * field[ll + (2 - 1), jj] * g0
374 field = temp_field
375
376 if time != -9999:
377Index: ../trunk-jpl/src/m/solve/parseresultsfromdisk.m
378===================================================================
379--- ../trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 25302)
380+++ ../trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 25303)
381@@ -228,15 +228,15 @@
382 g0 = md.love.g0;
383 mu0 = md.love.mu0;
384 rr = md.materials.radius;
385- rho= md.materials.density;
386- rho_avg = sum( rho .* diff(rr.^3)/sum(diff(rr.^3)) );
387- temp_field = cell(degmax+1,nfreq,nlayer+1,6);
388+ rho= md.materials.density;
389+ rho_avg = sum(rho.*diff(rr.^3)/sum(diff(rr.^3)));
390+ temp_field = cell(degmax+1,nfreq,nlayer+1,6);
391 for ii=1:degmax+1
392 for jj=1:nfreq
393 for kk=1:nlayer+1
394 if (kk<nlayer+1)
395- ll = (ii-1)*(nlayer+1)*6 + ((kk-1)*6+1) + 3;
396- temp_field{ii,jj,kk,1} = field(ll+(1-1),jj)*r0; % mm = 4
397+ ll = (ii-1)*(nlayer+1)*6 + ((kk-1)*6+1) + 3;
398+ temp_field{ii,jj,kk,1} = field(ll+(1-1),jj)*r0; % mm = 4
399 temp_field{ii,jj,kk,2} = field(ll+(2-1),jj)*mu0; % mm = 5
400 temp_field{ii,jj,kk,3} = field(ll+(3-1),jj)*r0; % mm = 6
401 temp_field{ii,jj,kk,4} = field(ll+(4-1),jj)*mu0; % mm = 1
402@@ -251,7 +251,7 @@
403 temp_field{ii,jj,kk,4} = 0;
404 if (md.love.forcing_type==9)
405 temp_field{ii,jj,kk,2} = 0;
406- temp_field{ii,jj,kk,6} = (2*ii-1)/r0 - ii*field(ll+(3-1),jj)*g0;
407+ temp_field{ii,jj,kk,6} = (2*ii-1)/r0 - ii*field(ll+(3-1),jj)*g0;
408 elseif (md.love.forcing_type==11)
409 temp_field{ii,jj,kk,2} = -(2*(ii-1)+1)*rho_avg/3;
410 temp_field{ii,jj,kk,6} = (2*ii-1)/r0 - ii*field(ll+(3-1),jj)*g0;
Note: See TracBrowser for help on using the repository browser.