Changeset 25129

Timestamp:

06/23/20 21:34:00 (5 years ago)

Author:

jdquinn

Message:

CHG: Saving progress on MATLAB to Python translations (still need to track down marshalling bug on Python side)

Location:

issm/trunk-jpl

Files:

• src/m/classes/giaivins.m (modified) (2 diffs)
• src/m/classes/giaivins.py (modified) (4 diffs)
• src/m/classes/giamme.m (modified) (2 diffs)
• src/m/classes/mesh3dsurface.m (modified) (1 diff)
• src/m/classes/mesh3dsurface.py (modified) (2 diffs)
• src/m/classes/model.py (modified) (2 diffs)
• src/m/parameterization/parameterize.py (modified) (1 diff)
• test/NightlyRun/test2001.py (modified) (2 diffs)
• test/NightlyRun/test2051.py (modified) (2 diffs)
• test/NightlyRun/test2052.py (modified) (1 diff)
• test/Par/GiaIvinsBenchmarksAB.py (modified) (1 diff)
• test/Par/GiaIvinsBenchmarksCD.py (modified) (1 diff)
• test/Par/SquareSheetConstrained.py (modified) (1 diff)

issm/trunk-jpl/src/m/classes/giaivins.m ¶

@@ -11 +11 @@
         end
         methods
-                function self = extrude(self,md) % {{{
-                        self.mantle_viscosity=project3d(md,'vector',self.mantle_viscosity,'type','node');
-                        self.lithosphere_thickness=project3d(md,'vector',self.lithosphere_thickness,'type','node');
-                end % }}}
                 function self = giaivins(varargin) % {{{
                         switch nargin
@@ -70 +66 @@
 
                 end % }}}
+                function self = extrude(self,md) % {{{
+                        self.mantle_viscosity=project3d(md,'vector',self.mantle_viscosity,'type','node');
+                        self.lithosphere_thickness=project3d(md,'vector',self.lithosphere_thickness,'type','node');
+                end % }}}
         end
 end

issm/trunk-jpl/src/m/classes/giaivins.py ¶

@@ +1 @@
+import numpy as np
+
+from checkfield import checkfield
 from fielddisplay import fielddisplay
 from project3d import project3d
-from checkfield import checkfield
 from WriteData import WriteData
 
@@ -7 +9 @@
 class giaivins(object):
     """
-    GIA class definition
+    GIA class definition for Ivins and James model
 
        Usage:
@@ -13 +15 @@
     """
 
-    def __init__(self):  # {{{
-        self.mantle_viscosity = float('NaN')
-        self.lithosphere_thickness = float('NaN')
-        self.cross_section_shape = 0
+    def __init__(self, *args): #{{{
+        self.mantle_viscosity       = np.nan
+        self.lithosphere_thickness  = np.nan
+        self.cross_section_shape    = 0
 
-    #set defaults
-        self.setdefaultparameters()
+        nargin = len(args)
 
+        if nargin == 0:
+            self.setdefaultparameters()
+        else:
+            raise Exception('constructor not supported')
     #}}}
 
-    def __repr__(self):  # {{{
-        string = '   giaivins solution parameters:'
+    def __repr__(self): #{{{
+        s = '   giaivins solution parameters:'
+        s += "{}\n".format(fielddisplay(self, 'mantle_viscosity', 'mantle viscosity constraints (NaN means no constraint) (Pa s)'))
+        s += "{}\n".format(fielddisplay(self, 'lithosphere_thickness', 'lithosphere thickness constraints (NaN means no constraint) (m)'))
+        s += "{}\n".format(fielddisplay(self, 'cross_section_shape', "1: square-edged, 2: elliptical-sedged surface"))
 
-        string = "%s\n%s" % (string, fielddisplay(self, 'mantle_viscosity', 'mantle viscosity constraints (NaN means no constraint) (Pa s)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'lithosphere_thickness', 'lithosphere thickness constraints (NaN means no constraint) (m)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'cross_section_shape', "1: square-edged, 2: elliptical - edged surface"))
-        return string
+        return s
     #}}}
 
-    def extrude(self, md):  # {{{
-        self.mantle_viscosity = project3d(md, 'vector', self.mantle_viscosity, 'type', 'node')
-        self.lithosphere_thickness = project3d(md, 'vector', self.lithosphere_thickness, 'type', 'node')
+    def setdefaultparameters(self): #{{{
+        self.cross_section_shape = 1
+
         return self
     #}}}
 
-    def setdefaultparameters(self):  # {{{
-        self.cross_section_shape = 1
-        return self
-    #}}}
-
-    def checkconsistency(self, md, solution, analyses):  # {{{
-        # Early return
-        if ('GiaAnalysis' not in analyses):
+    def checkconsistency(self, md, solution, analyses): #{{{
+        if not 'GiaAnalysis' in analyses:
             return md
 
@@ -56 +55 @@
 
         return md
-    # }}}
+    #}}}
 
-    def marshall(self, prefix, md, fid):  # {{{
-
+    def marshall(self, prefix, md, fid): #{{{
         WriteData(fid, prefix, 'object', self, 'fieldname', 'mantle_viscosity', 'format', 'DoubleMat', 'mattype', 1)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'lithosphere_thickness', 'format', 'DoubleMat', 'mattype', 1, 'scale', 10.**3.)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'cross_section_shape', 'format', 'Integer')
-    # }}}
+    #}}}
+
+    def extrude(self, md): #{{{
+        self.mantle_viscosity = project3d(md, 'vector', self.mantle_viscosity, 'type', 'node')
+        self.lithosphere_thickness = project3d(md, 'vector', self.lithosphere_thickness, 'type', 'node')
+        return self
+    #}}}

issm/trunk-jpl/src/m/classes/giamme.m ¶

@@ -14 +14 @@
         end
         methods
-                function self = extrude(self,md) % {{{
-                        for i=1:length(self.Ngia),
-                                self.Ngia{i}=project3d(md,'vector',self.Ngia{i},'type','node','layer',1);
-                                self.Ugia{i}=project3d(md,'vector',self.Ugia{i},'type','node','layer',1);
-                        end
-                end % }}}
                 function self = giamme(varargin) % {{{
                         switch nargin
@@ -74 +68 @@
 
                 end % }}}
+                function self = extrude(self,md) % {{{
+                        for i=1:length(self.Ngia),
+                                self.Ngia{i}=project3d(md,'vector',self.Ngia{i},'type','node','layer',1);
+                                self.Ugia{i}=project3d(md,'vector',self.Ugia{i},'type','node','layer',1);
+                        end
+                end % }}}
         end
 end

issm/trunk-jpl/src/m/classes/mesh3dsurface.m ¶

@@ -73 +73 @@
                 function obj = setdefaultparameters(obj) % {{{
 
-                        %the connectivity is the averaged number of nodes linked to a
-                        %given node through an edge. This connectivity is used to initially
-                        %allocate memory to the stiffness matrix. A value of 16 seems to
-                        %give a good memory/time ration. This value can be checked in
-                        %trunk/test/Miscellaneous/runme.m
+                        %The connectivity is the average number of nodes linked to a given 
+                        %node through an edge. This connectivity is used to initially 
+                        %allocate memory to the stiffness matrix. A value of 16 seems to 
+                        %give a good memory/time ratio. This value can be checked in 
+                        %test/NightlyRun/runme.py.
                         obj.average_vertex_connectivity=25;
                 end % }}}

issm/trunk-jpl/src/m/classes/mesh3dsurface.py ¶

@@ -9 +9 @@
 
 class mesh3dsurface(object):
-    #MESH3DSURFACE class definition
-    #
-    #   Usage:
-    #      mesh3dsurface = mesh3dsurface();
+    '''
+    MESH3DSURFACE class definition
+
+        Usage:
+            mesh3dsurface = mesh3dsurface()
+    '''
+
     def __init__(self, *args): # {{{
         self.x = np.nan
@@ -84 +87 @@
     # }}}
 
-    def loadobj(self): # {{{
-        # This def is directly called by matlab when a model() selfect is
-        # loaded. Update old properties here
-
-        #2014 Oct. 1st
-        if isstruct(self):
-            oldself = self
-            #Assign property values from struct
-            self = structtoobj(mesh3dsurface(), oldself)
-            if isfield(oldself, 'hemisphere'):
-                print('md.mesh.hemisphere has been automatically converted to EPSG code')
-                if strcmpi(oldself.hemisphere, 'n'):
-                    self.epsg = 3413
-                else:
-                    self.epsg = 3031
-        return self
-    # }}}
-
     def setdefaultparameters(self): # {{{
-        #the connectivity is the averaged number of nodes linked to a
-        #given node through an edge. This connectivity is used to initially
-        #allocate memory to the stiffness matrix. A value of 16 seems to
-        #give a good memory / time ration. This value can be checked in
-        #trunk / test / Miscellaneous / runme.m
+        #The connectivity is the average number of nodes linked to a given node 
+        #through an edge. This connectivity is used to initially allocate 
+        #memory to the stiffness matrix. A value of 16 seems to give a good 
+        #memory/time ratio. This value can be checked in 
+        #test/NightlyRun/runme.py.
         self.average_vertex_connectivity = 25
+
         return self
     # }}}

issm/trunk-jpl/src/m/classes/model.py ¶

@@ -54 +54 @@
 from transient import transient
 from giaivins import giaivins
+from giamme import giamme
 from esa import esa
 from autodiff import autodiff
@@ -117 +118 @@
         self.calving = calving()
         self.frontalforcings = frontalforcings()
-        self.gia = giaivins()
+        self.gia = giamme()
         self.love = fourierlove()
         self.esa = esa()

issm/trunk-jpl/src/m/parameterization/parameterize.py ¶

@@ -4 +4 @@
 
 def parameterize(md, parametername):
-    """
+    '''
     PARAMETERIZE - parameterize a model
 
-       from a parameter python file, start filling in all the model fields that were not
-       filled in by the mesh.py and mask.py model methods.
-       Warning: the parameter file must be able to be run in Python
+        From a parameter Python file, start filling in all the model fields 
+        that were not filled in by the mesh.py and mask.py model methods.
+        Warning: the parameter file must be able to be run in Python
 
-       Usage:
-          md = parameterize(md, parametername)
+        Usage:
+            md = parameterize(md, parametername)
 
-       Example:
-          md = parameterize(md, 'Square.par')
-    """
+        Example:
+            md = parameterize(md, 'Square.par')
+    '''
 
     #some checks

issm/trunk-jpl/test/NightlyRun/test2001.py ¶

@@ -4 +4 @@
 import numpy as np
 
+from giaivins import giaivins
 from model import *
 from parameterize import *
@@ -39 +40 @@
 #Solve for GIA deflection
 md.cluster = generic('name', gethostname(), 'np', 3)
+md.verbose = verbose('1111111')
 md = solve(md, 'Gia')
 

issm/trunk-jpl/test/NightlyRun/test2051.py ¶

@@ -1 +1 @@
 #Test Name: GiaIvinsBenchmarksAB2dA1
 import numpy as np
+
 from model import *
 from socket import gethostname
@@ -33 +34 @@
 field_names = ['GiaW', 'GiadWdt']
 field_tolerances = [1e-13, 1e-13]
-field_values = [md.results.GiaSolution.GiaW, md.results.GiaSolution.GiadWdt]
+field_values = [md.results.GiaSolution.UGia, md.results.GiaSolution.UGiaRate]

issm/trunk-jpl/test/NightlyRun/test2052.py ¶

@@ -1 +1 @@
 #Test Name: GiaIvinsBenchmarksAB2dA2
 import numpy as np
+
 from model import *
 from socket import gethostname

issm/trunk-jpl/test/Par/GiaIvinsBenchmarksAB.py ¶

@@ -1 +1 @@
 #Geometry specific to Experiments A and B
-
 import os.path
 import inspect
+
 import numpy as np
+
 from arch import *
+from giaivins import giaivins
 from InterpFromMeshToMesh2d import *
 from paterson import *

issm/trunk-jpl/test/Par/GiaIvinsBenchmarksCD.py ¶

@@ -31 +31 @@
 md.gia.mantle_viscosity = 1e21 * np.ones((md.mesh.numberofvertices, ))  #in Pa.s
 md.gia.lithosphere_thickness = 100 * np.ones((md.mesh.numberofvertices, ))  #in km
-md.materials.lithosphere_shear_modulus = 6.7 * 1e10  #in Pa
+md.materials.lithosphere_shear_modulus = 6.7e10  #in Pa
 md.materials.lithosphere_density = 3.32  #in g / cm^3
-md.materials.mantle_shear_modulus = 1.45 * 1e11  #in Pa
+md.materials.mantle_shear_modulus = 1.45e11  #in Pa
 md.materials.mantle_density = 3.34  #in g / cm^3
 

issm/trunk-jpl/test/Par/SquareSheetConstrained.py ¶

@@ -63 +63 @@
 #GIA:
 md.gia.lithosphere_thickness = 100. * np.ones((md.mesh.numberofvertices))  # in km
-md.gia.mantle_viscosity = 1. * 10**21 * np.ones((md.mesh.numberofvertices))  # in Pa.s
-md.materials.lithosphere_shear_modulus = 6.7 * 10**10  # in Pa
+md.gia.mantle_viscosity = 1.e21 * np.ones((md.mesh.numberofvertices))  # in Pa.s
+md.materials.lithosphere_shear_modulus = 6.7e10  # in Pa
 md.materials.lithosphere_density = 3.32  # in g / cm^ - 3
-md.materials.mantle_shear_modulus = 1.45 * 10**11  # in Pa
+md.materials.mantle_shear_modulus = 1.45e11  # in Pa
 md.materials.mantle_density = 3.34  # in g / cm^ - 3