Changeset 23431

Timestamp:

10/17/18 03:28:24 (6 years ago)

Author:

bdef

Message:

BUG: fix for model colapse

File:

• issm/trunk-jpl/src/m/classes/model.py (modified) (15 diffs)

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

@@ -16 +16 @@
 from basalforcings import basalforcings
 from matice import matice
-from levelset import levelset 
+from levelset import levelset
 from calving import calving
 from fourierlove import fourierlove
@@ -75 +75 @@
 
                 # classtype=model.properties
-                                
+
                 # for classe in dict.keys(classtype):
                 #       print classe
@@ -225 +225 @@
                    This routine extracts a submodel from a bigger model with respect to a given contour
                    md must be followed by the corresponding exp file or flags list
-                   It can either be a domain file (argus type, .exp extension), or an array of element flags. 
-                   If user wants every element outside the domain to be 
+                   It can either be a domain file (argus type, .exp extension), or an array of element flags.
+                   If user wants every element outside the domain to be
                    extract2d, add '~' to the name of the domain file (ex: '~HO.exp')
                    an empty string '' will be considered as an empty domain
@@ -349 +349 @@
                         md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos].astype(int)-1]
 
-                #Initial 2d mesh 
+                #Initial 2d mesh
                 if md1.mesh.__class__.__name__=='mesh3dprisms':
                         flag_elem_2d=flag_elem[np.arange(0,md1.mesh.numberofelements2d)]
@@ -431 +431 @@
                 if np.size(md1.stressbalance.spcvx)>1 and np.size(md1.stressbalance.spcvy)>2 and np.size(md1.stressbalance.spcvz)>2:
                         if np.size(md1.inversion.vx_obs)>1 and np.size(md1.inversion.vy_obs)>1:
-                                md2.stressbalance.spcvx[nodestoflag2]=md2.inversion.vx_obs[nodestoflag2] 
+                                md2.stressbalance.spcvx[nodestoflag2]=md2.inversion.vx_obs[nodestoflag2]
                                 md2.stressbalance.spcvy[nodestoflag2]=md2.inversion.vy_obs[nodestoflag2]
                         else:
@@ -508 +508 @@
                     - follow two polynomial laws, one for the lower part and one for the upper part of the mesh
                     - be discribed by a list of coefficients (between 0 and 1)
- 
+
 
                    Usage:
@@ -596 +596 @@
                 number_nodes3d=np.size(x3d)    #number of 3d nodes for the non extruded part of the mesh
 
-                #Extrude elements 
+                #Extrude elements
                 elements3d=np.empty((0,6),int)
                 for i in xrange(numlayers-1):
@@ -618 +618 @@
                 md.mesh.upperelements=upperelements
 
-                #Save old mesh 
+                #Save old mesh
                 md.mesh.x2d=md.mesh.x
                 md.mesh.y2d=md.mesh.y
@@ -625 +625 @@
                 md.mesh.numberofvertices2d=md.mesh.numberofvertices
 
-                #Build global 3d mesh 
+                #Build global 3d mesh
                 md.mesh.elements=elements3d
                 md.mesh.x=x3d
@@ -672 +672 @@
 
                 md.materials.extrude(md)
-                md.damage.extrude(md)
+                if md.damage.isdamage==1:
+                        md.damage.extrude(md)
                 md.gia.extrude(md)
                 md.mask.extrude(md)
@@ -688 +689 @@
                 '''
                 collapses a 3d mesh into a 2d mesh
-                        
+
                 This routine collapses a 3d model into a 2d model and collapses all
                 the fileds of the 3d model by taking their depth-averaged values
-                        
+
                 Usage:
                         md=collapse(md)
-                '''     
+                '''
 
                 #Check that the model is really a 3d model
                 if md.mesh.domaintype().lower() != '3d':
                         raise StandardError("only a 3D model can be collapsed")
-                
+
                 #dealing with the friction law
                 #drag is limited to nodes that are on the bedrock.
-                if hasattr(md.friction,'coefficient'): 
+                if hasattr(md.friction,'coefficient'):
                         md.friction.coefficient=project2d(md,md.friction.coefficient,1)
 
                 #p and q (same deal, except for element that are on the bedrock: )
-                if hasattr(md.friction,'p'): 
+                if hasattr(md.friction,'p'):
                         md.friction.p=project2d(md,md.friction.p,1)
-                if hasattr(md.friction,'q'): 
+                if hasattr(md.friction,'q'):
                         md.friction.q=project2d(md,md.friction.q,1)
-                
-                if hasattr(md.friction,'coefficientcoulomb'): 
+
+                if hasattr(md.friction,'coefficientcoulomb'):
                         md.friction.coefficientcoulomb=project2d(md,md.friction.coefficientcoulomb,1)
-                if hasattr(md.friction,'C'): 
+                if hasattr(md.friction,'C'):
                         md.friction.C=project2d(md,md.friction.C,1)
-                if hasattr(md.friction,'As'): 
+                if hasattr(md.friction,'As'):
                         md.friction.As=project2d(md,md.friction.As,1)
                 if hasattr(md.friction,'effective_pressure') and not np.isnan(md.friction.effective_pressure).all():
                         md.friction.effective_pressure=project2d(md,md.friction.effective_pressure,1)
-                if hasattr(md.friction,'water_layer'): 
+                if hasattr(md.friction,'water_layer'):
                         md.friction.water_layer=project2d(md,md.friction.water_layer,1)
-                if hasattr(md.friction,'m'): 
+                if hasattr(md.friction,'m'):
                         md.friction.m=project2d(md,md.friction.m,1)
 
                 #observations
-                if not np.isnan(md.inversion.vx_obs).all(): 
-                        md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers) 
-                if not np.isnan(md.inversion.vy_obs).all(): 
-                        md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers) 
-                if not np.isnan(md.inversion.vel_obs).all(): 
-                        md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers) 
-                if not np.isnan(md.inversion.cost_functions_coefficients).all(): 
-                        md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers) 
+                if not np.isnan(md.inversion.vx_obs).all():
+                        md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers)
+                if not np.isnan(md.inversion.vy_obs).all():
+                        md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers)
+                if not np.isnan(md.inversion.vel_obs).all():
+                        md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers)
+                if not np.isnan(md.inversion.cost_functions_coefficients).all():
+                        md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers)
                 if isinstance(md.inversion.min_parameters,np.ndarray):
-                        if md.inversion.min_parameters.size>1: 
-                                md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers) 
+                        if md.inversion.min_parameters.size>1:
+                                md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers)
                 if isinstance(md.inversion.max_parameters,np.ndarray):
-                        if md.inversion.max_parameters.size>1: 
-                                md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers) 
+                        if md.inversion.max_parameters.size>1:
+                                md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers)
                 if not np.isnan(md.smb.mass_balance).all():
-                        md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers) 
-                
+                        md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers)
+
                 #results
-                if not np.isnan(md.initialization.vx).all(): 
+                if not np.isnan(md.initialization.vx).all():
                         md.initialization.vx=DepthAverage(md,md.initialization.vx)
-                if not np.isnan(md.initialization.vy).all(): 
+                if not np.isnan(md.initialization.vy).all():
                         md.initialization.vy=DepthAverage(md,md.initialization.vy)
-                if not np.isnan(md.initialization.vz).all(): 
+                if not np.isnan(md.initialization.vz).all():
                         md.initialization.vz=DepthAverage(md,md.initialization.vz)
-                if not np.isnan(md.initialization.vel).all(): 
+                if not np.isnan(md.initialization.vel).all():
                         md.initialization.vel=DepthAverage(md,md.initialization.vel)
-                if not np.isnan(md.initialization.temperature).all(): 
+                if not np.isnan(md.initialization.temperature).all():
                         md.initialization.temperature=DepthAverage(md,md.initialization.temperature)
-                if not np.isnan(md.initialization.pressure).all(): 
+                if not np.isnan(md.initialization.pressure).all():
                         md.initialization.pressure=project2d(md,md.initialization.pressure,1)
-                if not np.isnan(md.initialization.sediment_head).all(): 
+                if not np.isnan(md.initialization.sediment_head).all():
                         md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1)
-                if not np.isnan(md.initialization.epl_head).all(): 
+                if not np.isnan(md.initialization.epl_head).all():
                         md.initialization.epl_head=project2d(md,md.initialization.epl_head,1)
-                if not np.isnan(md.initialization.epl_thickness).all(): 
+                if not np.isnan(md.initialization.epl_thickness).all():
                         md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1)
 
                 #giaivins
-                if not np.isnan(md.gia.mantle_viscosity).all(): 
-                        md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1) 
-                if not np.isnan(md.gia.lithosphere_thickness).all(): 
-                        md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1) 
+                if not np.isnan(md.gia.mantle_viscosity).all():
+                        md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1)
+                if not np.isnan(md.gia.lithosphere_thickness).all():
+                        md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1)
 
                 #elementstype
@@ -777 +778 @@
 
                 # Hydrologydc variables
-                if hasattr(md.hydrology,'hydrologydc'):
+                if type(md.hydrology is 'hydrologydc'):
                         md.hydrology.spcsediment_head=project2d(md,md.hydrology.spcsediment_head,1)
-                        md.hydrology.mask_eplactive_node=project2d(md,md.hydrology.mask_eplactive_node,1)
                         md.hydrology.sediment_transmitivity=project2d(md,md.hydrology.sediment_transmitivity,1)
                         md.hydrology.basal_moulin_input=project2d(md,md.hydrology.basal_moulin_input,1)
+                        md.hydrology.mask_thawed_node=project2d(md,md.hydrology.mask_thawed_node,1)
                         if md.hydrology.isefficientlayer == 1:
+                                md.hydrology.mask_eplactive_node=project2d(md,md.hydrology.mask_eplactive_node,1)
                                 md.hydrology.spcepl_head=project2d(md,md.hydrology.spcepl_head,1)
 
@@ -793 +795 @@
                 md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers)
                 md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers-1)
-                if not np.isnan(md.damage.spcdamage).all(): 
+                if not np.isnan(md.damage.spcdamage).all():
                         md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers-1)
 
@@ -799 +801 @@
                 md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B)
                 md.materials.rheology_n=project2d(md,md.materials.rheology_n,1)
-                
-                #damage: 
+
+                #damage:
                 if md.damage.isdamage:
                         md.damage.D=DepthAverage(md,md.damage.D)
 
                 #special for thermal modeling:
-                md.basalforcings.groundedice_melting_rate=project2d(md,md.basalforcings.groundedice_melting_rate,1) 
-                md.basalforcings.floatingice_melting_rate=project2d(md,md.basalforcings.floatingice_melting_rate,1) 
+                md.basalforcings.groundedice_melting_rate=project2d(md,md.basalforcings.groundedice_melting_rate,1)
+                md.basalforcings.floatingice_melting_rate=project2d(md,md.basalforcings.floatingice_melting_rate,1)
                 md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1) #bedrock only gets geothermal flux
 
@@ -848 +850 @@
                 mesh.numberofelements=md.mesh.numberofelements2d
                 mesh.elements=md.mesh.elements2d
-                if not np.isnan(md.mesh.vertexonboundary).all(): 
+                if not np.isnan(md.mesh.vertexonboundary).all():
                         mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1)
-                if not np.isnan(md.mesh.elementconnectivity).all(): 
+                if not np.isnan(md.mesh.elementconnectivity).all():
                         mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1)
                 if isinstance(md.mesh.lat,np.ndarray):
-                        if md.mesh.lat.size==md.mesh.numberofvertices:  
-                                mesh.lat=project2d(md,md.mesh.lat,1) 
+                        if md.mesh.lat.size==md.mesh.numberofvertices:
+                                mesh.lat=project2d(md,md.mesh.lat,1)
                 if isinstance(md.mesh.long,np.ndarray):
-                        if md.mesh.long.size==md.mesh.numberofvertices: 
-                                md.mesh.long=project2d(md,md.mesh.long,1) 
+                        if md.mesh.long.size==md.mesh.numberofvertices:
+                                md.mesh.long=project2d(md,md.mesh.long,1)
                 mesh.epsg=md.mesh.epsg
                 if isinstance(md.mesh.scale_factor,np.ndarray):
-                        if md.mesh.scale_factor.size==md.mesh.numberofvertices: 
-                                md.mesh.scale_factor=project2d(md,md.mesh.scale_factor,1) 
+                        if md.mesh.scale_factor.size==md.mesh.numberofvertices:
+                                md.mesh.scale_factor=project2d(md,md.mesh.scale_factor,1)
                 md.mesh=mesh
                 md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices)[0]