Ignore:
Timestamp:
06/01/22 05:01:48 (3 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 27033

Location:
issm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/test

  • issm/trunk/test/NightlyRun/test2084.m

    r26744 r27035  
    77
    88md=model();
    9 md.cluster=generic('name',oshostname(),'np',1);
     9md.cluster=generic('name',oshostname(),'np',8);
    1010
    11 % set validation=1 for comparing against the Spada benchark.
     11% set validation=1 for comparing against the Spada benchmark.
    1212validation=0;
    1313
    … …  
    1717
    1818md.verbose=verbose('all');
     19md.verbose=verbose('1111111111111111');
    1920cst=365.25*24*3600*1000;
    2021
    … …  
    3637
    3738md.love.allow_layer_deletion=1;
    38 md.love.frequencies=([0]*2*pi)'/cst;
     39md.love.frequencies=[0; (logspace(-6,3, 1000))'/cst];
    3940md.love.nfreq=length(md.love.frequencies);
    4041md.love.sh_nmin=1;
    41 md.love.sh_nmax=256;
     42md.love.sh_nmax=1000;
    4243md.love.underflow_tol=1e-20;
     44md.love.pw_threshold=1e-3;
    4345md.love.Gravitational_Constant=6.6732e-11;
    44 md.love.integration_steps_per_layer=200;
     46md.love.integration_steps_per_layer=100;
     47md.love.allow_layer_deletion=1;
     48md.love.forcing_type=11;
     49md.love.chandler_wobble=0;
     50md.love.complex_computation=0;
    4551
    4652md.love.istemporal=1;
    4753md.love.n_temporal_iterations=8;
    48 %md.love.time=(logspace(-4,5, 2))'*cst;
    49 md.love.time=(logspace(-1,2, 50))'*cst;
     54%md.love.time=(logspace(-6,5, 2))'*cst;
     55md.love.time=[0; (logspace(-3,5, 24))'*cst];
     56
     57%md.love.time=(linspace(1/12,10, 10*12))'*cst/1e3;
    5058md.love.love_kernels=1;
    5159if md.love.istemporal
    … …  
    5563md=solve(md,'lv');
    5664
    57 ht2=md.results.LoveSolution.LoveHr;
    58 lt2=md.results.LoveSolution.LoveLr;
    59 kt2=md.results.LoveSolution.LoveKr;
     65ht=md.results.LoveSolution.LoveHt';
     66lt=md.results.LoveSolution.LoveLt';
     67kt=md.results.LoveSolution.LoveKt';
    6068t=md.love.time/cst*1e3;
    6169
    … …  
    6977field_tolerances={2.0e-8,2.0e-8,2.0e-8};
    7078field_values={...
    71         (md.results.LoveSolution.LoveHr(:,1)),...
    72         (md.results.LoveSolution.LoveKr(:,1)),...
    73         (md.results.LoveSolution.LoveLr(:,1)),...
     79        (md.results.LoveSolution.LoveHt(:,1)),...
     80        (md.results.LoveSolution.LoveKt(:,1)),...
     81        (md.results.LoveSolution.LoveLt(:,1)),...
    7482        };
     83
     84return
     85
     86if false
     87addpath('../../../../invGIA/Spada_benchmark/')
     88load spada.mat
     89s_weak=[9 12 15];
     90
     91hspada(:,s_weak)=[];
     92kspada(:,s_weak)=[];
     93lspada(:,s_weak)=[];
     94
     95hts=zeros(length(t),md.love.sh_nmax+1);
     96lts=zeros(length(t),md.love.sh_nmax+1);
     97kts=zeros(length(t),md.love.sh_nmax+1);
     98for d=max(2,md.love.sh_nmin):md.love.sh_nmax
     99hts(:,d+1)=hspada(d-1,2);
     100lts(:,d+1)=lspada(d-1,2);
     101kts(:,d+1)=kspada(d-1,2);
     102for mo=1:9
     103hts(:,d+1)=hts(:,d+1)-hspada(d-1,3+mo)./sspada(d-1,1+mo).*(1-exp(t/1e3*sspada(d-1,1+mo)));
     104lts(:,d+1)=lts(:,d+1)-lspada(d-1,3+mo)./sspada(d-1,1+mo).*(1-exp(t/1e3*sspada(d-1,1+mo)));
     105kts(:,d+1)=kts(:,d+1)-kspada(d-1,3+mo)./sspada(d-1,1+mo).*(1-exp(t/1e3*sspada(d-1,1+mo)));
     106end
     107end
     108close all
     109
     110if md.love.sh_nmin==md.love.sh_nmax
     111subplot(2,3,1)
     112plot(t/cst*1e3,ht(:,md.love.sh_nmin+1),'x-',t/cst*1e3,hts);
     113set(gca, 'xscale', 'log'); shading flat;
     114title('h')
     115
     116subplot(2,3,2)
     117plot(t/cst*1e3,lt(:,md.love.sh_nmin+1),'x-',t/cst*1e3,lts);
     118set(gca, 'xscale', 'log'); shading flat;
     119title('l')
     120
     121subplot(2,3,3)
     122plot(t/cst*1e3,kt(:,md.love.sh_nmin+1),'x-',t/cst*1e3,kts);
     123set(gca, 'xscale', 'log'); shading flat;
     124title('k')
     125
     126subplot(2,3,4)
     127plot(t/cst*1e3,log10(abs((ht-hts)./hts))');
     128set(gca, 'xscale', 'log'); shading flat;
     129
     130subplot(2,3,5)
     131plot(t/cst*1e3,log10(abs((lt-lts)./lts))');
     132set(gca, 'xscale', 'log'); shading flat;
     133title('l')
     134
     135subplot(2,3,6)
     136plot(t/cst*1e3,log10(abs((kt-kts)./kts))');
     137set(gca, 'xscale', 'log'); shading flat;
     138title('k')
     139else
     140
     141[T,N]=meshgrid(t/cst*1e3,0:md.love.sh_nmax);
     142subplot(1,3,1)
     143pcolor(T,N,log10(abs((ht-hts)./hts))');
     144set(gca, 'xscale', 'log'); shading flat;colorbar;
     145title('h')
     146
     147subplot(1,3,2)
     148pcolor(T,N,log10(abs((lt-lts)./lts))');
     149set(gca, 'xscale', 'log'); shading flat;colorbar;
     150title('l')
     151
     152subplot(1,3,3)
     153pcolor(T,N,log10(abs((kt-kts)./kts))');
     154set(gca, 'xscale', 'log'); shading flat;colorbar;
     155title('k')
     156end
     157
    75158
    76159% validate elastic loading solutions against the Spada benchmark. {{{
    … …  
    100183        %
    101184end
     185
     186end
    102187% }}}
    103 
    104 md.love.frequencies=([1e-3 1e-2 1e-1 1 -1e-3 -1e-2 -1e-1 -1]*2*pi)'/cst;
    105 md.love.nfreq=length(md.love.frequencies);
    106 md.love.sh_nmax=256;
    107 md.materials.burgers_mu=md.materials.lame_mu;
    108 md.materials.burgers_viscosity=md.materials.viscosity;
    109 
    110 md=solve(md,'lv');
    111 
    112 %Fields and tolerances to track changes
    113 field_names     ={field_names{:},'LoveH_loading_realpart','LoveK_loading_realpart','LoveL_loading_realpart','LoveH_loading_imagpart','LoveK_loading_imagpart','LoveL_loading_imagpart'};
    114 field_tolerances={field_tolerances{:},5e-7,5e-7,5e-7,5e-7,5e-7,5e-7};
    115 field_values={field_values{:},...
    116         (md.results.LoveSolution.LoveHr(:,:)),...
    117         (md.results.LoveSolution.LoveKr(:,:)),...
    118         (md.results.LoveSolution.LoveLr(:,:)),...
    119         (md.results.LoveSolution.LoveHi(:,:)),...
    120         (md.results.LoveSolution.LoveKi(:,:)),...
    121         (md.results.LoveSolution.LoveLi(:,:)),...
    122         };
    123 
    124 md.love.forcing_type=9;
    125 md.love.sh_nmin=2;
    126 md.love.frequencies=([0 1e-3 1e-2 1e-1 1 -1e-3 -1e-2 -1e-1 -1]*2*pi)'/cst;
    127 md.love.nfreq=length(md.love.frequencies);
    128 md=solve(md,'lv');
    129 
    130 % validate elastic tidal solutions against the Spada benchmark. {{{
    131 if validation
    132         spada_solutions = load('spada_elastic_tidal_deg_h_l_k');
    133         spada_d = spada_solutions(:,1);
    134         spada_h = spada_solutions(:,2);
    135         spada_l = spada_solutions(:,3);
    136         spada_k = spada_solutions(:,4);
    137 
    138         %rename ISSM solutions. 
    139         issm_d = [2:md.love.sh_nmax]; 
    140         issm_h = md.results.LoveSolution.LoveHr(3:end,1);
    141         issm_l = md.results.LoveSolution.LoveLr(3:end,1);
    142         issm_k = md.results.LoveSolution.LoveKr(3:end,1);
    143 
    144         % relative difference for each degree, except for zero and one.
    145         diff_h = 1 - issm_h./spada_h; 
    146         diff_l = 1 - issm_l./spada_l; 
    147         diff_k = 1 - issm_k./spada_k; 
    148 
    149         % k seems to have larger relative difference at higher degrees, although values approach zero.
    150         figure
    151         plot(spada_d,issm_k./spada_k);
    152 
    153         figure
    154         plot(spada_d,[diff_h diff_l diff_k]); grid on;
    155         legend('h','l','k'); title('tidal love');
    156 else
    157         %
    158 end
    159 % }}}
    160 
    161 %tidal love numbers
    162 field_names     ={field_names{:},'LoveH_tidal_elastic','LoveK_tidal_elastic','LoveL_tidal_elastic','LoveH_tidal_realpart','LoveK_tidal_realpart','LoveL_tidal_realpart','LoveH_tidal_imagpart','LoveK_tidal_imagpart','LoveL_tidal_imagpart'};
    163 field_tolerances={field_tolerances{:},8e-6,8e-6,8e-6,8e-6,8e-6,8e-6,8e-6,8e-6,8e-6};
    164 field_values={field_values{:},...
    165         (md.results.LoveSolution.LoveHr(:,1)),...
    166         (md.results.LoveSolution.LoveKr(:,1)),...
    167         (md.results.LoveSolution.LoveLr(:,1)),...
    168         (md.results.LoveSolution.LoveHr(:,2:end)),...
    169         (md.results.LoveSolution.LoveKr(:,2:end)),...
    170         (md.results.LoveSolution.LoveLr(:,2:end)),...
    171         (md.results.LoveSolution.LoveHi(:,2:end)),...
    172         (md.results.LoveSolution.LoveKi(:,2:end)),...
    173         (md.results.LoveSolution.LoveLi(:,2:end)),...
    174         };
    175 
    176 %Many layers PREM-based model
    177 %data=load('../Data/PREM_500layers');
    178 %md.love.sh_nmin=1;
    179 %md.materials.radius=data(2:end-2,1);
    180 %md.materials.density=data(3:end-2,2);
    181 %md.materials.lame_lambda=data(3:end-2,3);
    182 %md.materials.lame_mu=data(3:end-2,4);
    183 %md.materials.issolid=data(3:end-2,4)>0;
    184 %ind=find(md.materials.issolid==0);
    185 %md.materials.density(ind(1))=sum((md.materials.radius(ind+1).^3-md.materials.radius(ind).^3).*md.materials.density(ind))/(md.materials.radius(ind(end)+1).^3-md.materials.radius(ind(1)+1).^3);
    186 %md.materials.lame_lambda(ind(1))=sum((md.materials.radius(ind+1).^3-md.materials.radius(ind).^3).*md.materials.lame_lambda(ind))/(md.materials.radius(ind(end)+1).^3-md.materials.radius(ind(1)+1).^3);
    187 %md.materials.lame_mu(ind(1))=sum((md.materials.radius(ind+1).^3-md.materials.radius(ind).^3).*md.materials.lame_mu(ind))/(md.materials.radius(ind(end)+1).^3-md.materials.radius(ind(1)).^3);
    188 %md.materials.radius(ind(2:end)+1)=[];
    189 %md.materials.density(ind(2:end))=[];
    190 %md.materials.lame_lambda(ind(2:end))=[];
    191 %md.materials.lame_mu(ind(2:end))=[];
    192 %md.materials.issolid(ind(2:end))=[];
    193 %md.materials.viscosity=10.^interp1([0 3479e3 3480e3 3680e3 5720e3 5800e3 6270e3 6371e3], log10([1e8 1e8 5e21 1e23 1e22 1e20 1e21 1e40]), md.materials.radius(2:end),'PCHIP');
    194 %md.materials.viscosity=md.materials.viscosity.*md.materials.issolid;
    195 %md.materials.burgers_mu=md.materials.lame_mu;
    196 %md.materials.burgers_viscosity=md.materials.viscosity;
    197 %md.materials.rheologymodel=md.materials.issolid*0;
    198 %md.love.forcing_type=11;
    199 %md.materials.numlayers=length(md.materials.viscosity);
    200 %md=solve(md,'lv');
    201 %
    202 %field_names     ={field_names{:},'LoveH_loadingVSS96_elastic','LoveK_loadingVSS96_elastic','LoveL_loadingVSS96_elastic','LoveH_loadingVSS96_realpart','LoveK_loadingVSS96_realpart','LoveL_loadingVSS96_realpart','LoveH_loadingVSS96_imagpart','LoveK_loadingVSS96_imagpart','LoveL_loadingVSS96_imagpart'};
    203 %field_tolerances={field_tolerances{:},2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6};
    204 %field_values={field_values{:},...
    205 %       (md.results.LoveSolution.LoveHr(:,1)),...
    206 %       (md.results.LoveSolution.LoveKr(:,1)),...
    207 %       (md.results.LoveSolution.LoveLr(:,1)),...
    208 %       (md.results.LoveSolution.LoveHr(:,2:end)),...
    209 %       (md.results.LoveSolution.LoveKr(:,2:end)),...
    210 %       (md.results.LoveSolution.LoveLr(:,2:end)),...
    211 %       (md.results.LoveSolution.LoveHi(:,2:end)),...
    212 %       (md.results.LoveSolution.LoveKi(:,2:end)),...
    213 %       (md.results.LoveSolution.LoveLi(:,2:end)),...
    214 %       };
    215 
    216 %Model VSS96 from Vermeersen, L.L.A., Sabadini, R. & Spada, G., 1996a. Analytical visco-elastic relaxation models, Geophys. Res. Lett., 23, 697–700.
    217 md.materials.radius=[10, 1222.5, 3480., 3600., 3630.5, 3700., 3900., 4000., 4200., 4300., 4500., 4600., 4800., 4900., 5100., 5200., 5400., 5500., 5600.5, 5650., 5701., 5736., 5771.5, 5821., 5951., 5970.5, 6016., 6061., 6150.5, 6151.5, 6251., 6371.]'*1e3;
    218 md.materials.lame_mu=[1e-5, 0., 2.933, 2.8990002, 2.8550003, 2.7340002, 2.675, 2.559, 2.502, 2.388, 2.331, 2.215, 2.157, 2.039, 1.979, 1.8560001, 1.794, 1.73, 1.639, 1.2390001, 1.224, 1.21, 1.128, 0.97700006, 0.906, 0.79, 0.773, 0.741, 0.656, 0.665, 0.602]'*1e11;
    219 md.materials.density=[10925., 10925., 5506.42, 5491.45, 5456.57, 5357.06, 5307.24, 5207.13, 5156.69, 5054.69, 5002.99, 4897.83, 4844.22, 4734.6, 4678.44, 4563.07, 4503.72, 4443.16, 4412.41, 3992.14, 3983.99, 3975.84, 3912.82, 3786.78, 3723.78, 3516.39, 3489.51, 3435.78, 3359.5, 3367.1, 3184.3]';
    220 md.materials.viscosity=[0., 0., 7.999999999999999E+21, 8.5E+21, 8.999999999999999E+21, 3.E+22, 4.E+22, 5.0000000000000004E+22, 6.E+22, 5.0000000000000004E+22, 4.5E+22, 3.E+22, 2.5000000000000002E+22, 1.7999999999999998E+22, 1.3E+22, 7.999999999999999E+21, 6.999999999999999E+21, 6.5E+21, 6.E+21, 5.5E+21, 5.E+21, 4.4999999999999995E+21, 3.9999999999999995E+21, 2.5E+21, 1.9999999999999997E+21, 1.5E+21, 9.999999999999999E+20, 6.E+20, 5.5000000000000007E+20, 2.E+20, 1.E40]';
    221 md.materials.lame_lambda=md.materials.lame_mu*0+5e14;
    222 md.materials.issolid=md.materials.lame_mu>0;
    223 md.materials.numlayers=length(md.materials.lame_mu);
    224 md.materials.burgers_mu=md.materials.lame_mu;
    225 md.materials.burgers_viscosity=md.materials.viscosity;
    226 md.materials.rheologymodel=md.materials.issolid*0;
    227 md.love.forcing_type=11;
    228 md.love.sh_nmin=1;
    229 md.love.sh_nmax=100;
    230 md=solve(md,'lv');
    231 md.love.frequencies=([0 1e-3 1e-2 1 -1e-3 -1e-2 -1]*2*pi)'/cst;
    232 md.love.nfreq=length(md.love.frequencies);
    233 
    234 field_names     ={field_names{:},'LoveH_loadingVSS96_elastic','LoveK_loadingVSS96_elastic','LoveL_loadingVSS96_elastic','LoveH_loadingVSS96_realpart','LoveK_loadingVSS96_realpart','LoveL_loadingVSS96_realpart','LoveH_loadingVSS96_imagpart','LoveK_loadingVSS96_imagpart','LoveL_loadingVSS96_imagpart'};
    235 field_tolerances={field_tolerances{:},2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6};
    236 field_values={field_values{:},...
    237         (md.results.LoveSolution.LoveHr(:,1)),...
    238         (md.results.LoveSolution.LoveKr(:,1)),...
    239         (md.results.LoveSolution.LoveLr(:,1)),...
    240         (md.results.LoveSolution.LoveHr(:,2:end)),...
    241         (md.results.LoveSolution.LoveKr(:,2:end)),...
    242         (md.results.LoveSolution.LoveLr(:,2:end)),...
    243         (md.results.LoveSolution.LoveHi(:,2:end)),...
    244         (md.results.LoveSolution.LoveKi(:,2:end)),...
    245         (md.results.LoveSolution.LoveLi(:,2:end)),...
    246         };
Note: See TracChangeset for help on using the changeset viewer.