Index: /issm/trunk-jpl/src/m/classes/geometry.py =================================================================== --- /issm/trunk-jpl/src/m/classes/geometry.py (revision 22095) +++ /issm/trunk-jpl/src/m/classes/geometry.py (revision 22096) @@ -47,4 +47,6 @@ if (solution=='TransientSolution' and md.transient.isgia) or (solution=='GiaSolution'): md = checkfield(md,'fieldname','geometry.thickness','NaN',1,'Inf',1,'>=',0,'timeseries',1) + elif solution=='LoveSolution': + return else: md = checkfield(md,'fieldname','geometry.surface' ,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) Index: /issm/trunk-jpl/src/m/classes/mask.py =================================================================== --- /issm/trunk-jpl/src/m/classes/mask.py (revision 22095) +++ /issm/trunk-jpl/src/m/classes/mask.py (revision 22096) @@ -38,4 +38,6 @@ #}}} def checkconsistency(self,md,solution,analyses): # {{{ + if(solution=='LoveSolution'): + return md = checkfield(md,'fieldname','mask.ice_levelset' ,'size',[md.mesh.numberofvertices]) Index: /issm/trunk-jpl/src/m/classes/materials.py =================================================================== --- /issm/trunk-jpl/src/m/classes/materials.py (revision 22095) +++ /issm/trunk-jpl/src/m/classes/materials.py (revision 22096) @@ -7,5 +7,5 @@ def naturetointeger(strnat): #{{{ - intnat=np.array(len(strnat)) + intnat=np.zeros(len(strnat)) for i in range(len(intnat)): if strnat[i]=='damageice': @@ -181,8 +181,8 @@ self.viscosity=[1e21,1e40] #mantle and lithosphere viscosity (respectively) [Pa.s] - self.lame_mu=[1.45*1e11,6.7*10^10] # (Pa) #lithosphere and mantle shear modulus (respectively) [Pa] + self.lame_mu=[1.45*1e11,6.7*1e10] # (Pa) #lithosphere and mantle shear modulus (respectively) [Pa] self.lame_lambda=self.lame_mu # (Pa) #mantle and lithosphere lamba parameter (respectively) [Pa] - self.burgers_viscosity=[NaN,NaN] - self.burgers_mu=[NaN,NaN] + self.burgers_viscosity=[np.nan,np.nan] + self.burgers_mu=[np.nan,np.nan] self.isburgers=[0,0] self.density=[5.51*1e3,5.50*1e3] # (Pa) #mantle and lithosphere density [kg/m^3] @@ -228,6 +228,6 @@ for i in range(md.materials.numlayers-1): - if md.materials.issolid[i] and md.materials.issolid[i+1]: #if there are at least two consecutive indices that contain issolid = 0 - raise RuntimeError("%s%i%s"%('2 or more adjacent fluid layers detected starting at layer',i,'This is not supported yet. Consider merging them.')) + if (not md.materials.issolid[i]) and (not md.materials.issolid[i+1]): #if there are at least two consecutive indices that contain issolid = 0 + raise RuntimeError("%s%i%s"%('2 or more adjacent fluid layers detected starting at layer ',i,'. This is not supported yet. Consider merging them.')) else: @@ -237,5 +237,5 @@ # }}} def marshall(self,prefix,md,fid): # {{{ - + #1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum WriteData(fid,prefix,'name','md.materials.type','data',6,'format','Integer') Index: /issm/trunk-jpl/src/m/classes/mesh2d.py =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh2d.py (revision 22095) +++ /issm/trunk-jpl/src/m/classes/mesh2d.py (revision 22096) @@ -82,4 +82,6 @@ #}}} def checkconsistency(self,md,solution,analyses): # {{{ + if(solution=='LoveSolution'): + return md = checkfield(md,'fieldname','mesh.x','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) Index: /issm/trunk-jpl/test/NightlyRun/test2084.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2084.py (revision 22096) +++ /issm/trunk-jpl/test/NightlyRun/test2084.py (revision 22096) @@ -0,0 +1,209 @@ +#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 * +from generic import generic +from materials import * + +md=model() +md.cluster=generic('name',gethostname(),'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=np.array([10,1222.5,3.4800e+03,5.7010e+03,5.9510e+03, + 6.3010e+03,6.3710e+03]).reshape(-1,1)*1e3 +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) +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 +md.materials.viscosity=np.array([0,0,2.0000e+00,1.0000e+00,1.0000e+00, + 1.0000e+25]).reshape(-1,1)*1e21 +md.materials.lame_lambda=np.array(md.materials.lame_mu)*0+5e14 +md.materials.issolid=np.array([1,0,1,1,1,1]).reshape(-1,1) +md.materials.isburgers=np.zeros((md.materials.numlayers,1)) + +md.love.allow_layer_deletion=1 +md.love.frequencies=(np.array([0])*2*pi).reshape(-1,1)/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=[3.2e-9,3.2e-9,3.2e-9] +field_values=[ +np.array(md.results.LoveSolution.LoveHr)[:,0], +np.array(md.results.LoveSolution.LoveKr)[:,0], +np.array(md.results.LoveSolution.LoveLr)[:,0] +] + +md.love.frequencies=(np.array([1e-3,1e-2,1e-1,1,-1e-3,-1e-2,-1e-1, + -1])*2*pi).reshape(-1,1)/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+=['LoveHlr','LoveKlr','LoveLlr','LoveHli','LoveKli','LoveLli'] +field_tolerances+=[3.2e-9,3.2e-9,3.2e-9,3.2e-9,3.2e-9,3.2e-9] +field_values+=[ +np.array(md.results.LoveSolution.LoveHr), +np.array(md.results.LoveSolution.LoveKr), +np.array(md.results.LoveSolution.LoveLr), +np.array(md.results.LoveSolution.LoveHi), +np.array(md.results.LoveSolution.LoveKi), +np.array(md.results.LoveSolution.LoveLi) +] + +md.love.forcing_type=9 +md.love.sh_nmin=2 +md.love.frequencies=((np.array([0,1e-3,1e-2,1e-1,1,-1e-3, + -1e-2,-1e-1,-1])*2*pi).reshape(-1,1)/cst) +md.love.nfreq=len(md.love.frequencies) + +md=solve(md,'lv') + + + +#tidal love numbers, check +field_names+=['LoveHtr','LoveKtr','LoveLtr','LoveHti','LoveKti','LoveLti'] +field_tolerances+=[3.2e-9,3.2e-9,3.2e-9,3.2e-9,3.2e-9,3.2e-9,3.2e-9,3.2e-9,3.2e-9] +field_values+=[ +np.array(md.results.LoveSolution.LoveHr)[:,0], +np.array(md.results.LoveSolution.LoveKr)[:,0], +np.array(md.results.LoveSolution.LoveLr)[:,0], +np.array(md.results.LoveSolution.LoveHr)[:,1:], +np.array(md.results.LoveSolution.LoveKr)[:,1:], +np.array(md.results.LoveSolution.LoveLr)[:,1:], +np.array(md.results.LoveSolution.LoveHi)[:,1:], +np.array(md.results.LoveSolution.LoveKi)[:,1:], +np.array(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=np.array([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.]).reshape(-1,1)*1e3 +md.materials.lame_mu=np.array([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]).reshape(-1,1)*1e11 +md.materials.density=np.array([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]).reshape(-1,1) +md.materials.viscosity=np.array([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]).reshape(-1,1) +md.materials.lame_lambda=np.array(md.materials.lame_mu)*0+5e14 +md.materials.issolid=np.ones(len(md.materials.lame_mu)).reshape(-1,1) +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=(np.array([0,1e-3,1e-2,1,-1e-3,-1e-2, + -1])*2*pi).reshape(-1,1)/cst +md.love.nfreq=len(md.love.frequencies) + +field_names+=['LoveHmr','LoveKmr','LoveLmr','LoveHmi','LoveKmi','LoveLmi'] +field_tolerances+=[3.2e-9,3.2e-9,3.2e-9,3.2e-9,3.2e-9,3.2e-9,3.2e-9,3.2e-9,3.2e-9] +field_values+=[ +np.array(md.results.LoveSolution.LoveHr)[:,0], +np.array(md.results.LoveSolution.LoveKr)[:,0], +np.array(md.results.LoveSolution.LoveLr)[:,0], +np.array(md.results.LoveSolution.LoveHr)[:,1:], +np.array(md.results.LoveSolution.LoveKr)[:,1:], +np.array(md.results.LoveSolution.LoveLr)[:,1:], +np.array(md.results.LoveSolution.LoveHi)[:,1:], +np.array(md.results.LoveSolution.LoveKi)[:,1:], +np.array(md.results.LoveSolution.LoveLi)[:,1:] +]