Index: /issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.m =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.m (revision 25165) +++ /issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.m (revision 25166) @@ -10062,5 +10062,5 @@ if strcmpi(frame,'CM'), return; - elseif strcmpi(frame,'CF'), % from Blewitt, 2003, JGR + elseif strcmpi(frame,'CF'), % from Blewitt, 2003, JGR if strcmpi(type,'loadingverticaldisplacement'), series(2,1) = -0.269; Index: /issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.py =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.py (revision 25165) +++ /issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.py (revision 25166) @@ -10069,9 +10069,9 @@ elif frame == 'CF': #from Blewitt, 2003, JGR if type == 'loadingverticaldisplacement': - series[1, 0] = -0.269 + series[0] = -0.269 elif type == 'loadinggravitationalpotential': - series[1, 0] = 0.021 + series[0] = 0.021 elif type == 'loadinghorizontaldisplacement': - series[1, 0] = 0.134 + series[0] = 0.134 else: raise Exception("love_numbers error message: unknown reference frame: {}".format(frame)) Index: /issm/trunk-jpl/src/m/classes/lovenumbers.py =================================================================== --- /issm/trunk-jpl/src/m/classes/lovenumbers.py (revision 25165) +++ /issm/trunk-jpl/src/m/classes/lovenumbers.py (revision 25166) @@ -61,5 +61,5 @@ #check that love numbers are provided at the same level of accuracy: - if (size(self.h, 0) != size(self.k, 0) | size(self.h, 0) != size(self.l, 0)): + if (self.h.shape[0] != self.k.shape[0]) or (self.h.shape[0] != self.l.shape[0]): raise ValueError('lovenumbers error message: love numbers should be provided at the same level of accuracy') #}}} Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.py (revision 25165) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.py (revision 25166) @@ -67,9 +67,9 @@ return md - md = checkfield(md, 'fieldname', 'solidearth.settings.reltol', 'size', [1, 1]) - md = checkfield(md, 'fieldname', 'solidearth.settings.abstol', 'size', [1, 1]) - md = checkfield(md, 'fieldname', 'solidearth.settings.maxiter', 'size', [1, 1], '>=', 1) - md = checkfield(md, 'fieldname', 'solidearth.settings.runfrequency', 'size', [1, 1], '>=', 1) - md = checkfield(md, 'fieldname', 'solidearth.settings.degacc', 'size', [1, 1], '>=', 1e-10) + md = checkfield(md, 'fieldname', 'solidearth.settings.reltol', 'size', [1]) + md = checkfield(md, 'fieldname', 'solidearth.settings.abstol', 'size', [1]) + md = checkfield(md, 'fieldname', 'solidearth.settings.maxiter', 'size', [1], '>=', 1) + md = checkfield(md, 'fieldname', 'solidearth.settings.runfrequency', 'size', [1], '>=', 1) + md = checkfield(md, 'fieldname', 'solidearth.settings.degacc', 'size', [1], '>=', 1e-10) md = checkfield(md, 'fieldname', 'solidearth.settings.horiz', 'NaN', 1, 'Inf', 1, 'values', [0, 1]) Index: /issm/trunk-jpl/test/NightlyRun/test2002.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2002.m (revision 25165) +++ /issm/trunk-jpl/test/NightlyRun/test2002.m (revision 25166) @@ -5,6 +5,6 @@ md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %500 km resolution mesh -%parameterize slr solution: -%slr loading: {{{ +%parameterize solidearth solution: +%solidearth loading: {{{ md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1); md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1); @@ -70,20 +70,26 @@ %eustatic + rigid run: -md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=0;md.solidearth.settings.rotation=0; +md.solidearth.settings.rigid=1; +md.solidearth.settings.elastic=0; +md.solidearth.settings.rotation=0; md=solve(md,'Sealevelrise'); Srigid=md.results.SealevelriseSolution.Sealevel; %eustatic + rigid + elastic run: -md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=1;md.solidearth.settings.rotation=0; +md.solidearth.settings.rigid=1; +md.solidearth.settings.elastic=1; +md.solidearth.settings.rotation=0; md=solve(md,'Sealevelrise'); Selastic=md.results.SealevelriseSolution.Sealevel; %eustatic + rigid + elastic + rotation run: -md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=1; +md.solidearth.settings.rigid=1; +md.solidearth.settings.elastic=1; +md.solidearth.settings.rotation=1; md=solve(md,'Sealevelrise'); Srotation=md.results.SealevelriseSolution.Sealevel; %Fields and tolerances to track changes -field_names ={'Eustatic','Rigid','Elastic','Rotation'}; +field_names={'Eustatic','Rigid','Elastic','Rotation'}; field_tolerances={1e-13,1e-13,1e-13,1e-13}; field_values={Seustatic,Srigid,Selastic,Srotation}; Index: /issm/trunk-jpl/test/NightlyRun/test2002.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2002.py (revision 25165) +++ /issm/trunk-jpl/test/NightlyRun/test2002.py (revision 25166) @@ -1,3 +1,3 @@ -#Test Name: EarthSolidearth +#Test Name: EarthSlr import numpy as np @@ -14,13 +14,13 @@ #mesh earth: md = model() -md.mesh = gmshplanet('radius', 6.371012e3, 'resolution', 700.) #500 km resolution mesh +md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 700.) #700 km resolution mesh #parameterize solidearth solution: #solidearth loading: -md.solidearth.surfaceload.icethicknesschange = np.zeros(md.mesh.numberofelements) -md.solidearth.sealevel = np.zeros(md.mesh.numberofvertices) +md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, 1)) +md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, 1)) md.dsl.global_average_thermosteric_sea_level_change = np.zeros((1, 1)) -md.dsl.sea_surface_height_change_above_geoid = np.zeros(md.mesh.numberofvertices + 1) -md.dsl.sea_water_pressure_change_at_sea_floor = np.zeros(md.mesh.numberofvertices + 1) +md.dsl.sea_surface_height_change_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1)) +md.dsl.sea_water_pressure_change_at_sea_floor = np.zeros((md.mesh.numberofvertices + 1, 1)) #antarctica @@ -29,5 +29,4 @@ pos = np.where(late < -80)[0] md.solidearth.surfaceload.icethicknesschange[pos] = -100 - #greenland pos = np.where(np.logical_and.reduce((late > 70, late < 80, longe > -60, longe < -30)))[0] @@ -36,12 +35,13 @@ #elastic loading from love numbers: md.solidearth.lovenumbers = lovenumbers('maxdeg', 100) +#}}} -#mask: +#mask: {{{ mask = gmtmask(md.mesh.lat, md.mesh.long) -icemask = np.ones((md.mesh.numberofvertices)) +icemask = np.ones(md.mesh.numberofvertices) pos = np.where(mask == 0)[0] icemask[pos] = -1 -pos = np.where(np.sum(mask[md.mesh.elements.astype(int) - 1], axis=1) < 3)[0] -icemask[md.mesh.elements[pos, :].astype(int) - 1] = -1 +pos = np.where(np.sum(mask[md.mesh.elements - 1], axis=1) < 3) +icemask[md.mesh.elements[pos, :] - 1] = -1 md.mask.ice_levelset = icemask md.mask.ocean_levelset = -icemask @@ -51,5 +51,5 @@ md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices)) -#make sure that the elements that have loads are fully grounded: +#make sure that the elements that have loads are fully grounded pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1 @@ -58,9 +58,10 @@ #pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] # Do we need to do this twice? md.mask.ice_levelset[md.mesh.elements[pos, :] - 1] = 1 +# }}} md.solidearth.settings.ocean_area_scaling = 0 #geometry for the bed, arbitrary -md.geometry.bed = -np.ones(md.mesh.numberofvertices) +md.geometry.bed = -np.ones((md.mesh.numberofvertices, 1)) #materials Index: /issm/trunk-jpl/test/NightlyRun/test2003.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2003.m (revision 25165) +++ /issm/trunk-jpl/test/NightlyRun/test2003.m (revision 25166) @@ -76,5 +76,6 @@ %eustatic + rigid + elastic + rotation run: md.solidearth.settings.rigid=1; -md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=1; +md.solidearth.settings.elastic=1; +md.solidearth.settings.rotation=1; md.cluster=generic('name',oshostname(),'np',3); %md.verbose=verbose('111111111'); Index: /issm/trunk-jpl/test/NightlyRun/test2003.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2003.py (revision 25165) +++ /issm/trunk-jpl/test/NightlyRun/test2003.py (revision 25166) @@ -1,3 +1,3 @@ -#Test Name: EarthSolidearth_rotationalFeedback +#Test Name: EarthSlr_rotationalFeedback import numpy as np Index: /issm/trunk-jpl/test/NightlyRun/test2101.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2101.m (revision 25165) +++ /issm/trunk-jpl/test/NightlyRun/test2101.m (revision 25166) @@ -8,5 +8,6 @@ % define load md.esa.deltathickness=zeros(md.mesh.numberofelements,1); -pos=450; md.esa.deltathickness(pos)=-100; % this is the only "icy" element +pos=450; +md.esa.deltathickness(pos)=-100; % this is the only "icy" element %love numbers: @@ -14,29 +15,31 @@ %mask: {{{ - md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); +md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); - %make sure wherever there is an ice load, that the mask is set to ice: - md.mask.ice_levelset=ones(md.mesh.numberofvertices,1); - pos=find(md.esa.deltathickness); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; +%make sure wherever there is an ice load, that the mask is set to ice: +md.mask.ice_levelset=ones(md.mesh.numberofvertices,1); +pos=find(md.esa.deltathickness); +md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; - %is ice grounded? - md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1); - pos=find(md.mask.ice_levelset<=0); md.mask.ocean_levelset(pos)=1; +%is ice grounded? +md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1); +pos=find(md.mask.ice_levelset<=0); +md.mask.ocean_levelset(pos)=1; % }}} %geometry: {{{ - di=md.materials.rho_ice/md.materials.rho_water; - md.geometry.thickness=ones(md.mesh.numberofvertices,1); - md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1); - md.geometry.base=md.geometry.surface-md.geometry.thickness; - md.geometry.bed=md.geometry.base; +di=md.materials.rho_ice/md.materials.rho_water; +md.geometry.thickness=ones(md.mesh.numberofvertices,1); +md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1); +md.geometry.base=md.geometry.surface-md.geometry.thickness; +md.geometry.bed=md.geometry.base; % }}} %materials: {{{ - md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); - md.materials.rheology_B=paterson(md.initialization.temperature); - md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); +md.materials.rheology_B=paterson(md.initialization.temperature); +md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); % }}} %Miscellaneous: {{{ - md.miscellaneous.name='test2101'; +md.miscellaneous.name='test2101'; % }}}