Changeset 25166

Timestamp:

06/26/20 20:45:42 (5 years ago)

Author:

jdquinn

Message:

CHG: Saving MATLAB -> Python changes; still hunting for implementation bug

Location:

issm/trunk-jpl

Files:

• src/m/boundaryconditions/getlovenumbers.m (modified) (1 diff)
• src/m/boundaryconditions/getlovenumbers.py (modified) (1 diff)
• src/m/classes/lovenumbers.py (modified) (1 diff)
• src/m/classes/solidearthsettings.py (modified) (1 diff)
• test/NightlyRun/test2002.m (modified) (2 diffs)
• test/NightlyRun/test2002.py (modified) (6 diffs)
• test/NightlyRun/test2003.m (modified) (1 diff)
• test/NightlyRun/test2003.py (modified) (1 diff)
• test/NightlyRun/test2101.m (modified) (2 diffs)

issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.m

@@ -10062 +10062 @@
         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; 

issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.py

@@ -10069 +10069 @@
     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))

issm/trunk-jpl/src/m/classes/lovenumbers.py

@@ -61 +61 @@
 
         #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')
     #}}}

issm/trunk-jpl/src/m/classes/solidearthsettings.py

@@ -67 +67 @@
             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])
 

issm/trunk-jpl/test/NightlyRun/test2002.m

@@ -5 +5 @@
 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 +70 @@
 
 %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};

issm/trunk-jpl/test/NightlyRun/test2002.py

@@ -1 @@
-#Test Name: EarthSolidearth
+#Test Name: EarthSlr
 import numpy as np
 
@@ -14 +14 @@
 #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 +29 @@
 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 +35 @@
 #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 +51 @@
 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 +58 @@
 #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

issm/trunk-jpl/test/NightlyRun/test2003.m

@@ -76 +76 @@
 %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');

issm/trunk-jpl/test/NightlyRun/test2003.py

@@ -1 @@
-#Test Name: EarthSolidearth_rotationalFeedback
+#Test Name: EarthSlr_rotationalFeedback
 import numpy as np
 

issm/trunk-jpl/test/NightlyRun/test2101.m

@@ -8 +8 @@
 % 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 +15 @@
 
 %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';
 % }}}