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

  • ../trunk-jpl/test/NightlyRun/test2085.m

     
    1010        md.materials.numlayers = 10;
    1111        md.love.forcing_type = 9;
    1212
    13         md.materials.density=  (1:md.materials.numlayers)'*0+5511;
    14         md.materials.lame_mu=  (1:md.materials.numlayers)'*0+0.75e11;
    15         md.materials.viscosity=(1:md.materials.numlayers)'*0+1e21;
    16         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;
    1717        md.materials.issolid=ones(md.materials.numlayers,1);
    1818        md.materials.isburgers=zeros(md.materials.numlayers,1);
    1919        md.materials.burgers_mu=md.materials.lame_mu/3;
    … …  
    2424
    2525        md.love.allow_layer_deletion=1;
    2626        md.love.frequencies=0;
    27         md.love.nfreq=length(md.love.frequencies);
     27        md.love.nfreq=length(md.love.frequencies); % TODO: Why are we grabbing length of a scalar here?
    2828
    2929        md.love.sh_nmin = 2;
    3030        md.love.sh_nmax = 2;
    … …  
    3737        md=solve(md,'lv');
    3838
    3939        % extract love kernels {{{
    40         y1=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,1))); 
     40        y1=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,1)));
    4141        y2=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,2)));
    4242        y3=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,3)));
    4343        y4=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,4)));
    … …  
    147147        'y6_load_degree2_interface1','y6_load_degree2_interface2','y6_load_degree2_interface7','y6_load_degree2_interface10','y6_load_degree2_interface11',...
    148148        };
    149149field_tolerances={...
    150         1e-10,1e-10,1e-10,1e-10,1e-10,...
    151         1e-10,1e-10,1e-10,1e-10,1e-10,...
    152         1e-10,1e-10,1e-10,1e-10,1e-10,...
    153         1e-10,1e-10,1e-10,1e-10,1e-10,...
    154         1e-10,1e-10,1e-10,1e-10,1e-10,...
    155         1e-10,1e-10,1e-10,1e-10,1e-10,...
    156         1e-10,1e-10,1e-10,1e-10,1e-10,...
    157         1e-10,1e-10,1e-10,1e-10,1e-10,...
    158         1e-10,1e-10,1e-10,1e-10,1e-10,...
    159         1e-10,1e-10,1e-10,1e-10,1e-10,...
    160         1e-10,1e-10,1e-10,1e-10,1e-10,...
    161         1e-10,1e-10,1e-10,1e-10,1e-10,...
     150        9e-8, 1e-7, 3e-7, 3e-7, 3e-7,...
     151        9e-8, 1e-7, 3e-7, 3e-7, 1e-10,...
     152        9e-8, 8e-8, 2e-8, 2e-7, 4e-7,...
     153        9e-8, 9e-8, 2e-7, 4e-7, 1e-10,...
     154        4e-7, 4e-7, 2e-7, 3e-8, 2e-8,...
     155        2e-5, 2e-6, 2e-6, 1e-6, 2e-7,...
     156        3e-6, 3e-6, 3e-6, 4e-6, 4e-6,...
     157        3e-6, 3e-6, 2e-6, 6e-7, 1e-10,...
     158        3e-6, 3e-6, 5e-7, 3e-6, 5e-6,...
     159        3e-6, 3e-6, 9e-7, 7e-7, 1e-10,...
     160        4e-6, 4e-6, 3e-6, 5e-7, 3e-7,...
     161        3e-6, 3e-6, 2e-6, 7e-7, 2e-7...
    162162        };
    163163field_values={...
    164164        y1_tidal_degree2_interface1, y1_tidal_degree2_interface2, y1_tidal_degree2_interface7, y1_tidal_degree2_interface10,y1_tidal_degree2_interface11,...

  • ../trunk-jpl/test/NightlyRun/test2085.py

     
    1 #Test Name: LovenumberstAtDepth.
    2 #Same as test  #1 of test2084.m
     1# Test Name: FourierLoveKernels
     2# Homogenous Earth, for which analytic solutions exist.
     3# Love kernels for degree 2 (tested against analytic solns).
    34
    4 from model import *
     5
    56from socket import gethostname
    6 from solve import *
    7 from numpy import *
     7
     8import numpy as np
     9
    810from generic import generic
    911from materials import *
     12from model import *
     13from solve import *
    1014
     15# Fro volumetric potential
    1116md = model()
    12 md.cluster = generic('name', gethostname(), 'np', 1)
    1317
    1418md.materials = materials('litho')
    15 md.miscellaneous.name = 'FourierLoveTest'
    16 md.groundingline.migration = 'None'
    1719
    18 md.verbose = verbose('111111101')
    19 cst = 365.25 * 24 * 3600 * 1000
     20md.materials.numlayers = 10
     21md.love.forcing_type = 9
    2022
    21 md.materials.numlayers = 6
    22 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
    23 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)
    24 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
    25 md.materials.viscosity = np.array([0, 0, 2.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+25]).reshape(-1, 1) * 1e21
    26 md.materials.lame_lambda = np.array(md.materials.lame_mu) * 0 + 5e14
    27 md.materials.issolid = np.array([1, 0, 1, 1, 1, 1]).reshape(-1, 1)
    28 md.materials.isburgers = np.zeros((md.materials.numlayers)).reshape(-1, 1)
     23md.materials.density = np.zeros((md.materials.numlayers, 1)) + 5511
     24md.materials.lame_mu = np.zeros((md.materials.numlayers, 1)) + 0.75e11
     25md.materials.viscosity = np.zeros((md.materials.numlayers, 1)) + 1e21
     26md.materials.lame_lambda = np.zeros((md.materials.numlayers, 1)) + 0.5e17
     27md.materials.issolid = np.ones((md.materials.numlayers, 1))
     28md.materials.isburgers = np.zeros((md.materials.numlayers, 1))
     29md.materials.burgers_mu = md.materials.lame_mu / 3
     30md.materials.burgers_viscosity = md.materials.viscosity / 10
    2931
    30 md.love.love_kernels = 1
     32md.materials.radius = np.linspace(10e3, 6371e3, md.materials.numlayers + 1).reshape(-1, 1)
     33md.love.g0 = 9.8134357285509388 # directly grabbed from fourierlovesolver for this particular case
     34
    3135md.love.allow_layer_deletion = 1
    32 md.love.frequencies = (np.array([0]) * 2 * pi).reshape(-1, 1) / cst
    33 md.love.nfreq = len(md.love.frequencies)
     36md.love.frequencies = 0
     37md.love.nfreq = 1
     38
     39md.love.sh_nmin = 2
    3440md.love.sh_nmax = 2
     41md.love.love_kernels = 1
    3542
    36 md.materials.burgers_mu = md.materials.lame_mu
    37 md.materials.burgers_viscosity = md.materials.viscosity
     43md.miscellaneous.name = 'kernels'
     44md.cluster = generic('name', gethostname(), 'np', 1)
     45md.verbose = verbose('111111101')
    3846
    39 md = solve(md, 'lv')
     47md = solve(md,'lv')
    4048
     49# Extract love kernels #{{{
     50y1 = md.results.LoveSolution.LoveKernelsReal[:,0,:,0].squeeze()
     51y2 = md.results.LoveSolution.LoveKernelsReal[:,0,:,1].squeeze()
     52y3 = md.results.LoveSolution.LoveKernelsReal[:,0,:,2].squeeze()
     53y4 = md.results.LoveSolution.LoveKernelsReal[:,0,:,3].squeeze()
     54y5 = md.results.LoveSolution.LoveKernelsReal[:,0,:,4].squeeze()
     55y6 = md.results.LoveSolution.LoveKernelsReal[:,0,:,5].squeeze()
     56
     57y1_tidal_degree2_interface1 = y1[2, 0]
     58y1_tidal_degree2_interface2 = y1[2, 1]
     59y1_tidal_degree2_interface7 = y1[2, 6]
     60y1_tidal_degree2_interface10 = y1[2, 9]
     61y1_tidal_degree2_interface11 = y1[2, 10]
     62
     63y2_tidal_degree2_interface1 = y2[2, 0]
     64y2_tidal_degree2_interface2 = y2[2, 1]
     65y2_tidal_degree2_interface7 = y2[2, 6]
     66y2_tidal_degree2_interface10 = y2[2, 9]
     67y2_tidal_degree2_interface11 = y2[2, 10]
     68
     69y3_tidal_degree2_interface1 = y3[2, 0]
     70y3_tidal_degree2_interface2 = y3[2, 1]
     71y3_tidal_degree2_interface7 = y3[2, 6]
     72y3_tidal_degree2_interface10 = y3[2, 9]
     73y3_tidal_degree2_interface11 = y3[2, 10]
     74
     75y4_tidal_degree2_interface1 = y4[2, 0]
     76y4_tidal_degree2_interface2 = y4[2, 1]
     77y4_tidal_degree2_interface7 = y4[2, 6]
     78y4_tidal_degree2_interface10 = y4[2, 9]
     79y4_tidal_degree2_interface11 = y4[2, 10]
     80
     81y5_tidal_degree2_interface1 = y5[2, 0]
     82y5_tidal_degree2_interface2 = y5[2, 1]
     83y5_tidal_degree2_interface7 = y5[2, 6]
     84y5_tidal_degree2_interface10 = y5[2, 9]
     85y5_tidal_degree2_interface11 = y5[2, 10]
     86
     87y6_tidal_degree2_interface1 = y6[2, 0]
     88y6_tidal_degree2_interface2 = y6[2, 1]
     89y6_tidal_degree2_interface7 = y6[2, 6]
     90y6_tidal_degree2_interface10 = y6[2, 9]
     91y6_tidal_degree2_interface11 = y6[2, 10]
     92#}}}
     93
     94# For surface load
     95md.love.forcing_type = 11
     96
     97md = solve(md,'lv')
     98
     99# Extract love kernels #{{{
     100y1 = md.results.LoveSolution.LoveKernelsReal[:,0,:,0].squeeze()
     101y2 = md.results.LoveSolution.LoveKernelsReal[:,0,:,1].squeeze()
     102y3 = md.results.LoveSolution.LoveKernelsReal[:,0,:,2].squeeze()
     103y4 = md.results.LoveSolution.LoveKernelsReal[:,0,:,3].squeeze()
     104y5 = md.results.LoveSolution.LoveKernelsReal[:,0,:,4].squeeze()
     105y6 = md.results.LoveSolution.LoveKernelsReal[:,0,:,5].squeeze()
     106
     107y1_load_degree2_interface1 = y1[2, 0]
     108y1_load_degree2_interface2 = y1[2, 1]
     109y1_load_degree2_interface7 = y1[2, 6]
     110y1_load_degree2_interface10 = y1[2, 9]
     111y1_load_degree2_interface11 = y1[2, 10]
     112
     113y2_load_degree2_interface1 = y2[2, 0]
     114y2_load_degree2_interface2 = y2[2, 1]
     115y2_load_degree2_interface7 = y2[2, 6]
     116y2_load_degree2_interface10 = y2[2, 9]
     117y2_load_degree2_interface11 = y2[2, 10]
     118
     119y3_load_degree2_interface1 = y3[2, 0]
     120y3_load_degree2_interface2 = y3[2, 1]
     121y3_load_degree2_interface7 = y3[2, 6]
     122y3_load_degree2_interface10 = y3[2, 9]
     123y3_load_degree2_interface11 = y3[2, 10]
     124
     125y4_load_degree2_interface1 = y4[2, 0]
     126y4_load_degree2_interface2 = y4[2, 1]
     127y4_load_degree2_interface7 = y4[2, 6]
     128y4_load_degree2_interface10 = y4[2, 9]
     129y4_load_degree2_interface11 = y4[2, 10]
     130
     131y5_load_degree2_interface1 = y5[2, 0]
     132y5_load_degree2_interface2 = y5[2, 1]
     133y5_load_degree2_interface7 = y5[2, 6]
     134y5_load_degree2_interface10 = y5[2, 9]
     135y5_load_degree2_interface11 = y5[2, 10]
     136
     137y6_load_degree2_interface1 = y6[2, 0]
     138y6_load_degree2_interface2 = y6[2, 1]
     139y6_load_degree2_interface7 = y6[2, 6]
     140y6_load_degree2_interface10 = y6[2, 9]
     141y6_load_degree2_interface11 = y6[2, 10]
     142#}}}
     143
     144
    41145#Fields and tolerances to track changes
    42146#loading love numbers
    43 field_names = ['LoveH_loading_elastic', 'LoveK_loading_elastic', 'LoveL_loading_elastic', 'LoveKernels_degree1', 'LoveKernels_degree2']
    44 field_tolerances = [1e-10, 1e-10, 1e-10, 1e-10, 1e-10]
    45 field_values = [md.results.LoveSolution.LoveHr[:, 0],
    46                 md.results.LoveSolution.LoveKr[:, 0],
    47                 md.results.LoveSolution.LoveLr[:, 0],
    48                 md.results.LoveSolution.LoveKernelsReal[1, 0, :, :],
    49                 md.results.LoveSolution.LoveKernelsReal[2, 0, :, :]]
     147field_names = [
     148    'y1_tidal_degree2_interface1', 'y1_tidal_degree2_interface2', 'y1_tidal_degree2_interface7', 'y1_tidal_degree2_interface10', 'y1_tidal_degree2_interface11',
     149    'y2_tidal_degree2_interface1', 'y2_tidal_degree2_interface2', 'y2_tidal_degree2_interface7', 'y2_tidal_degree2_interface10', 'y2_tidal_degree2_interface11',
     150    'y3_tidal_degree2_interface1', 'y3_tidal_degree2_interface2', 'y3_tidal_degree2_interface7', 'y3_tidal_degree2_interface10', 'y3_tidal_degree2_interface11',
     151    'y4_tidal_degree2_interface1', 'y4_tidal_degree2_interface2', 'y4_tidal_degree2_interface7', 'y4_tidal_degree2_interface10', 'y4_tidal_degree2_interface11',
     152    'y5_tidal_degree2_interface1', 'y5_tidal_degree2_interface2', 'y5_tidal_degree2_interface7', 'y5_tidal_degree2_interface10', 'y5_tidal_degree2_interface11',
     153    'y6_tidal_degree2_interface1', 'y6_tidal_degree2_interface2', 'y6_tidal_degree2_interface7', 'y6_tidal_degree2_interface10', 'y6_tidal_degree2_interface11',
     154    'y1_load_degree2_interface1', 'y1_load_degree2_interface2', 'y1_load_degree2_interface7', 'y1_load_degree2_interface10', 'y1_load_degree2_interface11',
     155    'y2_load_degree2_interface1', 'y2_load_degree2_interface2', 'y2_load_degree2_interface7', 'y2_load_degree2_interface10', 'y2_load_degree2_interface11',
     156    'y3_load_degree2_interface1', 'y3_load_degree2_interface2', 'y3_load_degree2_interface7', 'y3_load_degree2_interface10', 'y3_load_degree2_interface11',
     157    'y4_load_degree2_interface1', 'y4_load_degree2_interface2', 'y4_load_degree2_interface7', 'y4_load_degree2_interface10', 'y4_load_degree2_interface11',
     158    'y5_load_degree2_interface1', 'y5_load_degree2_interface2', 'y5_load_degree2_interface7', 'y5_load_degree2_interface10', 'y5_load_degree2_interface11',
     159    'y6_load_degree2_interface1', 'y6_load_degree2_interface2', 'y6_load_degree2_interface7', 'y6_load_degree2_interface10', 'y6_load_degree2_interface11'
     160    ]
     161field_tolerances = [
     162    9e-8, 1e-7, 3e-7, 3e-7, 3e-7,
     163    9e-8, 1e-7, 3e-7, 3e-7, 1e-10,
     164    9e-8, 8e-8, 2e-8, 2e-7, 4e-7,
     165    9e-8, 9e-8, 2e-7, 4e-7, 1e-10,
     166    4e-7, 4e-7, 2e-7, 3e-8, 2e-8,
     167    2e-5, 2e-6, 2e-6, 1e-6, 2e-7,
     168    3e-6, 3e-6, 3e-6, 4e-6, 4e-6,
     169    3e-6, 3e-6, 2e-6, 6e-7, 1e-10,
     170    3e-6, 3e-6, 5e-7, 3e-6, 5e-6,
     171    3e-6, 3e-6, 9e-7, 7e-7, 1e-10,
     172    4e-6, 4e-6, 3e-6, 5e-7, 3e-7,
     173    3e-6, 3e-6, 2e-6, 7e-7, 2e-7
     174    ]
     175field_values = [
     176    y1_tidal_degree2_interface1, y1_tidal_degree2_interface2, y1_tidal_degree2_interface7, y1_tidal_degree2_interface10,y1_tidal_degree2_interface11,
     177    y2_tidal_degree2_interface1, y2_tidal_degree2_interface2, y2_tidal_degree2_interface7, y2_tidal_degree2_interface10,y2_tidal_degree2_interface11,
     178    y3_tidal_degree2_interface1, y3_tidal_degree2_interface2, y3_tidal_degree2_interface7, y3_tidal_degree2_interface10,y3_tidal_degree2_interface11,
     179    y4_tidal_degree2_interface1, y4_tidal_degree2_interface2, y4_tidal_degree2_interface7, y4_tidal_degree2_interface10,y4_tidal_degree2_interface11,
     180    y5_tidal_degree2_interface1, y5_tidal_degree2_interface2, y5_tidal_degree2_interface7, y5_tidal_degree2_interface10,y5_tidal_degree2_interface11,
     181    y6_tidal_degree2_interface1, y6_tidal_degree2_interface2, y6_tidal_degree2_interface7, y6_tidal_degree2_interface10,y6_tidal_degree2_interface11,
     182    y1_load_degree2_interface1, y1_load_degree2_interface2, y1_load_degree2_interface7, y1_load_degree2_interface10,y1_load_degree2_interface11,
     183    y2_load_degree2_interface1, y2_load_degree2_interface2, y2_load_degree2_interface7, y2_load_degree2_interface10,y2_load_degree2_interface11,
     184    y3_load_degree2_interface1, y3_load_degree2_interface2, y3_load_degree2_interface7, y3_load_degree2_interface10,y3_load_degree2_interface11,
     185    y4_load_degree2_interface1, y4_load_degree2_interface2, y4_load_degree2_interface7, y4_load_degree2_interface10,y4_load_degree2_interface11,
     186    y5_load_degree2_interface1, y5_load_degree2_interface2, y5_load_degree2_interface7, y5_load_degree2_interface10,y5_load_degree2_interface11,
     187    y6_load_degree2_interface1, y6_load_degree2_interface2, y6_load_degree2_interface7, y6_load_degree2_interface10,y6_load_degree2_interface11
     188    ]

  • ../trunk-jpl/test/NightlyRun/test2052.m

     
    3030%% solve for GIA deflection
    3131md=solve(md,'Gia');
    3232
    33 pause
    34 
    3533%Test Name: GiaIvinsBenchmarksAB2dC1
    3634U_AB2dC1 = md.results.GiaSolution.UGia;
    3735URate_AB2dC1 = md.results.GiaSolution.UGiaRate;

  • ../trunk-jpl/src/m/classes/geometry.m

     
    8989
    9090                end % }}}
    9191                function marshall(self,prefix,md,fid) % {{{
    92                         disp(md.geometry.thickness)
    9392                        WriteData(fid,prefix,'object',self,'fieldname','surface','format','DoubleMat','mattype',1);
    9493                        WriteData(fid,prefix,'object',self,'fieldname','thickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    95                         disp(md.geometry.thickness)
    96                         pause
    9794                        WriteData(fid,prefix,'object',self,'fieldname','base','format','DoubleMat','mattype',1);
    9895                        WriteData(fid,prefix,'object',self,'fieldname','bed','format','DoubleMat','mattype',1);
    9996                        WriteData(fid,prefix,'object',self,'fieldname','hydrostatic_ratio','format','DoubleMat','mattype',1);

  • ../trunk-jpl/src/m/solve/parseresultsfromdisk.py

     
    247247            mu0 = md.love.mu0
    248248            rr = md.materials.radius
    249249            rho = md.materials.density
    250             rho_avg = (rho * np.diff(np.power(rr, 3), n=1, axis=0)) / np.diff(np.power(rr, 3).sum()).sum()
     250            rho_avg_partial = np.diff(np.power(rr, 3), n=1, axis=0)
     251            rho_avg = ((rho * rho_avg_partial) / rho_avg_partial.sum()).sum()
    251252            temp_field = np.empty((degmax + 1, nfreq, nlayer + 1, 6))
    252253            temp_field.fill(0.0)
    253254            for ii in range(degmax + 1):
    254255                for jj in range(nfreq):
    255256                    for kk in range(nlayer + 1):
    256                         if kk < nlayer + 1:
     257                        if kk < nlayer: # NOTE: Upper bound of range is non-inclusive (compare to src/m/solve/parseresultsfromdisk.m)
    257258                            ll = ii * (nlayer + 1) * 6 + (kk * 6 + 1) + 3
    258                             temp_field[ii, jj, kk, 0] = field[ll + (1 - 1), jj] * r0        # mm = 4
    259                             temp_field[ii, jj, kk, 1] = field[ll + (2 - 1), jj] * mu0       # mm = 5
    260                             temp_field[ii, jj, kk, 2] = field[ll + (3 - 1), jj] * r0        # mm = 6
    261                             temp_field[ii, jj, kk, 3] = field[ll + (4 - 1), jj] * mu0       # mm = 1
    262                             temp_field[ii, jj, kk, 4] = field[ll + (5 - 1), jj] * r0 * g0   # mm = 2
    263                             temp_field[ii, jj, kk, 5] = field[ll + (6 - 1), jj] * g0        # mm = 3
    264                             print(temp_field)
     259                            temp_field[ii, jj, kk, 0] = field[ll + (0 - 1), jj] * r0        # mm = 4
     260                            temp_field[ii, jj, kk, 1] = field[ll + (1 - 1), jj] * mu0       # mm = 5
     261                            temp_field[ii, jj, kk, 2] = field[ll + (2 - 1), jj] * r0        # mm = 6
     262                            temp_field[ii, jj, kk, 3] = field[ll + (3 - 1), jj] * mu0       # mm = 1
     263                            temp_field[ii, jj, kk, 4] = field[ll + (4 - 1), jj] * r0 * g0   # mm = 2
     264                            temp_field[ii, jj, kk, 5] = field[ll + (5 - 1), jj] * g0        # mm = 3
    265265                        else: # surface
    266                             ll = ii * (nlayer + 1) * 6 - 2
    267                             temp_field[ii, jj, kk, 0] = field[ll + (1 - 1), jj] * r0
    268                             temp_field[ii, jj, kk, 2] = field[ll + (2 - 1), jj] * r0
    269                             temp_field[ii, jj, kk, 4] = field[ll + (3 - 1), jj] * r0 * g0
     266                            ll = (ii + 1) * (nlayer + 1) * 6 - 2
     267                            temp_field[ii, jj, kk, 0] = field[ll + (0 - 1), jj] * r0
     268                            temp_field[ii, jj, kk, 2] = field[ll + (1 - 1), jj] * r0
     269                            temp_field[ii, jj, kk, 4] = field[ll + (2 - 1), jj] * r0 * g0
    270270                            # surface BC
    271271                            temp_field[ii, jj, kk, 3] = 0
    272272                            if md.love.forcing_type == 9:
    273273                                temp_field[ii, jj, kk, 1] = 0
    274                                 temp_field[ii, jj, kk, 5] = (2 * ii - 1) / r0 - ii * field[ll + (3 - 1), jj] * g0
     274                                temp_field[ii, jj, kk, 5] = (2 * (ii + 1) - 1) / r0 - (ii + 1) * field[ll + (2 - 1), jj] * g0
    275275                            elif md.love.forcing_type == 11:
    276                                 temp_field[ii, jj, kk, 1] = -(2 * (ii - 1) + 1) * rho_avg / 3
    277                                 temp_field[ii, jj, kk, 5] = (2 * ii - 1) / r0 - ii * field[ll + (3 - 1), jj] * g0
     276                                temp_field[ii, jj, kk, 1] = -(2 * ii + 1) * rho_avg / 3
     277                                temp_field[ii, jj, kk, 5] = (2 * (ii + 1) - 1) / r0 - (ii + 1) * field[ll + (2 - 1), jj] * g0
    278278            field = temp_field
    279279
    280280        if time != -9999:

  • ../trunk-jpl/src/m/solve/parseresultsfromdisk.m

     
    228228                g0 = md.love.g0;
    229229                mu0 = md.love.mu0;
    230230                rr = md.materials.radius;
    231                 rho= md.materials.density; 
    232                 rho_avg = sum( rho .*  diff(rr.^3)/sum(diff(rr.^3)) );
    233                 temp_field = cell(degmax+1,nfreq,nlayer+1,6); 
     231                rho= md.materials.density;
     232                rho_avg = sum(rho.*diff(rr.^3)/sum(diff(rr.^3)));
     233                temp_field = cell(degmax+1,nfreq,nlayer+1,6);
    234234                for ii=1:degmax+1
    235235                        for jj=1:nfreq
    236236                                for kk=1:nlayer+1
    237237                                        if (kk<nlayer+1)
    238                                                 ll = (ii-1)*(nlayer+1)*6 + ((kk-1)*6+1) + 3; 
    239                                                 temp_field{ii,jj,kk,1} = field(ll+(1-1),jj)*r0;         % mm = 4 
     238                                                ll = (ii-1)*(nlayer+1)*6 + ((kk-1)*6+1) + 3;
     239                                                temp_field{ii,jj,kk,1} = field(ll+(1-1),jj)*r0;         % mm = 4
    240240                                                temp_field{ii,jj,kk,2} = field(ll+(2-1),jj)*mu0;        % mm = 5
    241241                                                temp_field{ii,jj,kk,3} = field(ll+(3-1),jj)*r0;         % mm = 6
    242242                                                temp_field{ii,jj,kk,4} = field(ll+(4-1),jj)*mu0;        % mm = 1
    … …  
    251251                                                temp_field{ii,jj,kk,4} = 0;
    252252                                                if (md.love.forcing_type==9)
    253253                                                        temp_field{ii,jj,kk,2} = 0;
    254                                                         temp_field{ii,jj,kk,6} = (2*ii-1)/r0 - ii*field(ll+(3-1),jj)*g0; 
     254                                                        temp_field{ii,jj,kk,6} = (2*ii-1)/r0 - ii*field(ll+(3-1),jj)*g0;
    255255                                                elseif (md.love.forcing_type==11)
    256256                                                        temp_field{ii,jj,kk,2} = -(2*(ii-1)+1)*rho_avg/3;
    257257                                                        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.