Index: ../trunk-jpl/test/NightlyRun/test2084.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test2084.m (nonexistent) +++ ../trunk-jpl/test/NightlyRun/test2084.m (revision 22060) @@ -0,0 +1,162 @@ + +%Test Name: GiaCaron. +%Forward Love number solution for a viscoelastic earth, model M3-L70-V01 from Spada, G., Barletta, V. R., Klemann, V., Riva, R. E. M., Martinec, Z., Gasperini, P., Lund, B., Wolf, D., Vermeersen, L. L. A. and King, M. A. (2011), A benchmark study for glacial isostatic adjustment codes. Geophysical Journal International, 185: 106–132. doi:10.1111/j.1365-246X.2011.04952.x + +md=model(); +md.cluster=generic('name',oshostname(),'np',1); + +md.materials=materials('litho'); +md.miscellaneous.name='FourierLoveTest'; + +md.verbose=verbose('111111101'); +cst=365.25*24*3600*1000; + + md.materials.numlayers=6; + md.materials.radius = [10 1222.5 3.4800e+03 5.7010e+03 5.9510e+03 6.3010e+03 6.3710e+03]'*1e3; + md.materials.density= [1.0750e4 1.0750e+04 4.9780e+03 3.8710e+03 3.4380e+03 3.0370e+03]'; + md.materials.lame_mu= [1e-5 0 2.2834e+00 1.0549e+00 7.0363e-01 5.0605e-01]'*1e11; + md.materials.viscosity=[0 0 2.0000e+00 1.0000e+00 1.0000e+00 1.0000e+25]'*1e21; + md.materials.lame_lambda=md.materials.lame_mu*0+5e14; + md.materials.issolid=[1 0 1 1 1 1]'; + md.materials.isburgers=zeros(md.materials.numlayers,1); + +md.love.allow_layer_deletion=1; +md.love.frequencies=([0]*2*pi)'/cst; +md.love.nfreq=length(md.love.frequencies); +md.love.sh_nmax=256; + +md.materials.burgers_mu=md.materials.lame_mu; +md.materials.burgers_viscosity=md.materials.viscosity; + +md=solve(md,'lv'); + + + +%Fields and tolerances to track changes +%loading love numbers +field_names ={'LoveH_loading_elastic','LoveK_loading_elastic','LoveL_loading_elastic'}; +field_tolerances={1e-13,1e-13,1e-13}; +field_values={... + (md.results.LoveSolution.LoveHr(:,1)),... + (md.results.LoveSolution.LoveKr(:,1)),... + (md.results.LoveSolution.LoveLr(:,1)),... + }; + + +md.love.frequencies=([1e-3 1e-2 1e-1 1 -1e-3 -1e-2 -1e-1 -1]*2*pi)'/cst; +md.love.nfreq=length(md.love.frequencies); +md.love.sh_nmax=256; +md.materials.burgers_mu=md.materials.lame_mu; +md.materials.burgers_viscosity=md.materials.viscosity; + +md=solve(md,'lv'); +%Fields and tolerances to track changes +%loading love numbers +field_names ={field_names{:},'LoveH_loading_realpart','LoveK_loading_realpart','LoveL_loading_realpart','LoveH_loading_imagpart','LoveK_loading_imagpart','LoveL_loading_imagpart'}; +field_tolerances={field_tolerances{:},1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; +field_values={field_values{:},... + (md.results.LoveSolution.LoveHr(:,:)),... + (md.results.LoveSolution.LoveKr(:,:)),... + (md.results.LoveSolution.LoveLr(:,:)),... + (md.results.LoveSolution.LoveHi(:,:)),... + (md.results.LoveSolution.LoveKi(:,:)),... + (md.results.LoveSolution.LoveLi(:,:)),... + }; + + +return +md.love.forcing_type=9; +md.love.sh_nmin=2 +md.love.frequencies=([0 1e-4 1e-3 1e-2 1e-1 1 -1e-4 -1e-3 -1e-2 -1e-1 -1]*2*pi)'/cst; +md.love.nfreq=length(md.love.frequencies); +md=solve(md,'lv'); +%tidal love numbers +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'}; +field_tolerances={field_tolerances{:},1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; +field_values={field_values{:},... + (md.results.LoveSolution.LoveHr(:,1)),... + (md.results.LoveSolution.LoveKr(:,1)),... + (md.results.LoveSolution.LoveLr(:,1)),... + (md.results.LoveSolution.LoveHr(:,2:end)),... + (md.results.LoveSolution.LoveKr(:,2:end)),... + (md.results.LoveSolution.LoveLr(:,2:end)),... + (md.results.LoveSolution.LoveHi(:,2:end)),... + (md.results.LoveSolution.LoveKi(:,2:end)),... + (md.results.LoveSolution.LoveLi(:,2:end)),... + }; + + + + %Many layers PREM-based model +%data=load('../Data/PREM_500layers'); +%md.love.sh_nmin=1; +%md.materials.radius=data(2:end-2,1); +%md.materials.density=data(3:end-2,2); +%md.materials.lame_lambda=data(3:end-2,3); +%md.materials.lame_mu=data(3:end-2,4); +%md.materials.issolid=data(3:end-2,4)>0; +%ind=find(md.materials.issolid==0); +%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); +%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); +%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); +%md.materials.radius(ind(2:end)+1)=[]; +%md.materials.density(ind(2:end))=[]; +%md.materials.lame_lambda(ind(2:end))=[]; +%md.materials.lame_mu(ind(2:end))=[]; +%md.materials.issolid(ind(2:end))=[]; +%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'); +%md.materials.viscosity=md.materials.viscosity.*md.materials.issolid; +%md.materials.burgers_mu=md.materials.lame_mu; +%md.materials.burgers_viscosity=md.materials.viscosity; +%md.materials.isburgers=md.materials.issolid*0; +%md.love.forcing_type=11; +%md.materials.numlayers=length(md.materials.viscosity); +%md=solve(md,'lv'); +% +%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'}; +%field_tolerances={field_tolerances{:},1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; +%field_values={field_values{:},... +% (md.results.LoveSolution.LoveHr(:,1)),... +% (md.results.LoveSolution.LoveKr(:,1)),... +% (md.results.LoveSolution.LoveLr(:,1)),... +% (md.results.LoveSolution.LoveHr(:,2:end)),... +% (md.results.LoveSolution.LoveKr(:,2:end)),... +% (md.results.LoveSolution.LoveLr(:,2:end)),... +% (md.results.LoveSolution.LoveHi(:,2:end)),... +% (md.results.LoveSolution.LoveKi(:,2:end)),... +% (md.results.LoveSolution.LoveLi(:,2:end)),... +% }; + + +%Model VSS96 from Vermeersen, L.L.A., Sabadini, R. & Spada, G., 1996a. Analytical visco-elastic relaxation models, Geophys. Res. Lett., 23, 697–700. +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; +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; +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]'; +md.materials.viscosity=[0., 0., 8E+21, 8.5E+21, 9E+21, 3.E+22, 4.E+22, 5.E+22, 6.E+22, 5.E+22, 4.5E+22, 3.E+22, 2.5E+22, 1.8E+22, 1.3E+22, 8E+21, 7+21, 6.5E+21, 6.E+21, 5.5E+21, 5.E+21, 4.5E+21, 4E+21, 2.5E+21, 2.E+21, 1.5E+21, 10.E+20, 6.E+20, 5.5E+20, 2.E+20, 1.E40]'; +md.materials.lame_lambda=md.materials.lame_mu*0+5e14; +md.materials.issolid=md.materials.lame_mu>0; +md.materials.numlayers=length(md.materials.lame_mu); +md.materials.burgers_mu=md.materials.lame_mu; +md.materials.burgers_viscosity=md.materials.viscosity; +md.materials.isburgers=md.materials.issolid*0; +md.love.forcing_type=11; +md.love.sh_nmin=1; +md.love.sh_nmax=100; +md=solve(md,'lv'); +md.love.frequencies=([0 1e-4 1e-2 1 -1e-4 -1e-2 -1]*2*pi)'/cst; +md.love.nfreq=length(md.love.frequencies); + + +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'}; +field_tolerances={field_tolerances{:},1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; +field_values={field_values{:},... + (md.results.LoveSolution.LoveHr(:,1)),... + (md.results.LoveSolution.LoveKr(:,1)),... + (md.results.LoveSolution.LoveLr(:,1)),... + (md.results.LoveSolution.LoveHr(:,2:end)),... + (md.results.LoveSolution.LoveKr(:,2:end)),... + (md.results.LoveSolution.LoveLr(:,2:end)),... + (md.results.LoveSolution.LoveHi(:,2:end)),... + (md.results.LoveSolution.LoveKi(:,2:end)),... + (md.results.LoveSolution.LoveLi(:,2:end)),... + }; Index: ../trunk-jpl/test/NightlyRun/test2084.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test2084.py (nonexistent) +++ ../trunk-jpl/test/NightlyRun/test2084.py (revision 22060) @@ -0,0 +1,152 @@ + +#Test Name: GiaCaron. Forward Love number solution for a viscoelastic earth, model M3-L70-V01 from Spada, G., Barletta, V. R., Klemann, V., Riva, R. E. M., Martinec, Z., Gasperini, P., Lund, B., Wolf, D., Vermeersen, L. L. A. and King, M. A. (2011), A benchmark study for glacial isostatic adjustment codes. Geophysical Journal International, 185: 106–132. doi:10.1111/j.1365-246X.2011.04952.x +from model import * +from socket import gethostname +from solve import * +from numpy import * +md=model() +md.cluster=generic('name',oshostname(),'np',1) +md.materials=materials('litho') +md.miscellaneous.name='FourierLoveTest' +md.verbose=verbose('111111101') +cst=365.25*24*3600*1000 +md.materials.numlayers=6 +md.materials.radius =[10, 1222.5, 3.4800e+03,5.7010e+03,5.9510e+03,6.3010e+03,6.3710e+03]*1e3 +md.materials.density=[1.0750e4, 1.0750e+04,4.9780e+03,3.8710e+03,3.4380e+03,3.0370e+03] +md.materials.lame_mu=[1e-5,0,2.2834e+00,1.0549e+00,7.0363e-01,5.0605e-01]*1e11 +md.materials.viscosity=[0,0,2.0000e+00,1.0000e+00,1.0000e+00,1.0000e+25]*1e21 +md.materials.lame_lambda=md.materials.lame_mu*0+5e14 +md.materials.issolid=[1, 0, 1, 1, 1, 1] +md.materials.isburgers=zeros(md.materials.numlayers,1) +md.love.allow_layer_deletion=1 + +md.love.frequencies=([0]*2*pi)/cst +md.love.nfreq=len(md.love.frequencies) +md.love.sh_nmax=256 +md.materials.burgers_mu=md.materials.lame_mu +md.materials.burgers_viscosity=md.materials.viscosity +md=solve(md,'lv') +#Fields and tolerances to track changes +#loading love numbers +field_names=['LoveHe','LoveKe','LoveLe'] +field_tolerances=[1e-13,1e-13,1e-13] +field_values=[\ +(md.results.LoveSolution.LoveHr[:][0]),\ +(md.results.LoveSolution.LoveKr[:][0]),\ +(md.results.LoveSolution.LoveLr[:][0]),\ +] + + +md.love.frequencies=([1e-3, 1e-2, 1e-1, 1, -1e-3, -1e-2, -1e-1, -1]*2*pi)/cst +md.love.nfreq=len(md.love.frequencies) +md.love.sh_nmax=256 +md.materials.burgers_mu=md.materials.lame_mu +md.materials.burgers_viscosity=md.materials.viscosity +md=solve(md,'lv') +#Fields and tolerances to track changes +#loading love numbers +field_names=concatenate([field_names,['LoveHlr','LoveKlr','LoveLlr','LoveHli','LoveKli','LoveLli']]) +field_tolerances=concatenate([field_tolerances,[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13]]) +field_values=concatenate([field_values,[\ +(md.results.LoveSolution.LoveHr[:][:]),\ +(md.results.LoveSolution.LoveKr[:][:]),\ +(md.results.LoveSolution.LoveLr[:][:])),\ +(md.results.LoveSolution.LoveHi[:][:]),\ +(md.results.LoveSolution.LoveKi[:][:]),\ +(md.results.LoveSolution.LoveLi[:][:]),\ +]]) + + + +md.love.forcing_type=9 +md.love.sh_nmin=2 +md.love.frequencies=([0, 1e-4, 1e-3, 1e-2, 1e-1, 1, -1e-4, -1e-3, -1e-2, -1e-1, -1]*2*pi)/cst +md.love.nfreq=len(md.love.frequencies) +md=solve(md,'lv') +#tidal love numbers +field_names=concatenate([field_names,['LoveHtr','LoveKtr','LoveLtr','LoveHti','LoveKti','LoveLti']]) +field_tolerances=concatenate([field_tolerances,[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13]]) +field_values=concatenate([field_values,[\ +(md.results.LoveSolution.LoveHr[:][0]),\ +(md.results.LoveSolution.LoveKr[:][0]),\ +(md.results.LoveSolution.LoveLr[:][0]),\ +(md.results.LoveSolution.LoveHr[:][1:]),\ +(md.results.LoveSolution.LoveKr[:][1:]),\ +(md.results.LoveSolution.LoveLr[:][1:]),\ +(md.results.LoveSolution.LoveHi[:][1:]),\ +(md.results.LoveSolution.LoveKi[:][1:]),\ +(md.results.LoveSolution.LoveLi[:][1:]),\ +]]) + + + +#Many layers PREM-based model +#data=load('../Data/PREM_500layers') +#md.love.sh_nmin=1 +#md.materials.radius=data(2:end-2,1) +#md.materials.density=data(3:end-2,2) +#md.materials.lame_lambda=data(3:end-2,3) +#md.materials.lame_mu=data(3:end-2,4) +#md.materials.issolid=data(3:end-2,4)>0 +#ind=find(md.materials.issolid==0) +#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) +#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) +#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) +#md.materials.radius(ind(2:end)+1)=[] +#md.materials.density(ind(2:end))=[] +#md.materials.lame_lambda(ind(2:end))=[] +#md.materials.lame_mu(ind(2:end))=[] +#md.materials.issolid(ind(2:end))=[] +#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') +#md.materials.viscosity=md.materials.viscosity.*md.materials.issolid +#md.materials.burgers_mu=md.materials.lame_mu +#md.materials.burgers_viscosity=md.materials.viscosity +#md.materials.isburgers=md.materials.issolid*0 +#md.love.forcing_type=11 +#md.materials.numlayers=len(md.materials.viscosity) +#md=solve(md,'lv') +# +#field_names=[field_names,'LoveHmr','LoveKmr','LoveLmr','LoveHmi','LoveKmi','LoveLmi'] +#field_tolerances=[field_tolerances,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13] +#field_values=[field_values,\ +# (md.results.LoveSolution.LoveHr[:][0]),\ +# (md.results.LoveSolution.LoveKr[:][0]),\ +# (md.results.LoveSolution.LoveLr[:][0]),\ +# (md.results.LoveSolution.LoveHr[:][1:]),\ +# (md.results.LoveSolution.LoveKr[:][1:]),\ +# (md.results.LoveSolution.LoveLr[:][1:]),\ +# (md.results.LoveSolution.LoveHi[:][1:]),\ +# (md.results.LoveSolution.LoveKi[:][1:]),\ +# (md.results.LoveSolution.LoveLi[:][1:]),\ +# ] +#Model VSS96 from Vermeersen, L.L.A., Sabadini, R. & Spada, G., 1996a. Analytical visco-elastic relaxation models, Geophys. Res. Lett., 23, 697–700. +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 +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 +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] +md.materials.viscosity=[0., 0., 8E+21, 8.5E+21, 9E+21, 3.E+22, 4.E+22, 5E+22, 6.E+22, 5E+22, 4.5E+22, 3.E+22, 2.5E+22, 1.8E+22, 1.3E+22, 8E+21, 7E+21, 6.5E+21, 6.E+21, 5.5E+21, 5.E+21, 55E+21, 45E+21, 2.5E+21, 27E+21, 1.5E+21, 10E+20, 6.E+20, 5.5E+20, 2.E+20, 1.E40] +md.materials.lame_lambda=md.materials.lame_mu*0+5e14 +md.materials.issolid=ones(len(md.materials.mu)) +md.materials.issolid[1]=0 +md.materials.numlayers=len(md.materials.lame_mu) +md.materials.burgers_mu=md.materials.lame_mu +md.materials.burgers_viscosity=md.materials.viscosity +md.materials.isburgers=md.materials.issolid*0 +md.love.forcing_type=11 +md.love.sh_nmin=1 +md.love.sh_nmax=100 +md=solve(md,'lv') +md.love.frequencies=([0, 1e-4, 1e-2, 1, -1e-4, -1e-2, -1]*2*pi)/cst +md.love.nfreq=len(md.love.frequencies) +field_names=concatenate([field_names,['LoveHmr','LoveKmr','LoveLmr','LoveHmi','LoveKmi','LoveLmi']]) +field_tolerances=concatenate([field_tolerances,[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13]]) +field_values=concatenate([field_values,[\ +(md.results.LoveSolution.LoveHr[:][0]),\ +(md.results.LoveSolution.LoveKr[:][0]),\ +(md.results.LoveSolution.LoveLr[:][0]),\ +(md.results.LoveSolution.LoveHr[:][1:]),\ +(md.results.LoveSolution.LoveKr[:][1:]),\ +(md.results.LoveSolution.LoveLr[:][1:]),\ +(md.results.LoveSolution.LoveHi[:][1:]),\ +(md.results.LoveSolution.LoveKi[:][1:]),\ +(md.results.LoveSolution.LoveLi[:][1:]),\ +]])