Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 22056) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 22057) @@ -65,5 +65,4 @@ from ElementConnectivity import ElementConnectivity from contourenvelope import contourenvelope -import MatlabFuncs as m from DepthAverage import DepthAverage #}}} @@ -126,43 +125,43 @@ def properties(self): # {{{ # ordered list of properties since vars(self) is random - return ['mesh',\ - 'mask',\ - 'geometry',\ - 'constants',\ - 'smb',\ - 'basalforcings',\ - 'materials',\ - 'damage',\ - 'friction',\ - 'flowequation',\ - 'timestepping',\ - 'initialization',\ - 'rifts',\ - 'slr',\ - 'debug',\ - 'verbose',\ - 'settings',\ - 'toolkits',\ - 'cluster',\ - 'balancethickness',\ - 'stressbalance',\ - 'groundingline',\ - 'hydrology',\ - 'masstransport',\ - 'thermal',\ - 'steadystate',\ - 'transient',\ - 'levelset',\ - 'calving',\ - 'gia',\ - 'love',\ - 'autodiff',\ - 'inversion',\ - 'qmu',\ - 'amr',\ - 'outputdefinition',\ - 'results',\ - 'radaroverlay',\ - 'miscellaneous',\ + return ['mesh', + 'mask', + 'geometry', + 'constants', + 'smb', + 'basalforcings', + 'materials', + 'damage', + 'friction', + 'flowequation', + 'timestepping', + 'initialization', + 'rifts', + 'slr', + 'debug', + 'verbose', + 'settings', + 'toolkits', + 'cluster', + 'balancethickness', + 'stressbalance', + 'groundingline', + 'hydrology', + 'masstransport', + 'thermal', + 'steadystate', + 'transient', + 'levelset', + 'calving', + 'gia', + 'love', + 'autodiff', + 'inversion', + 'qmu', + 'amr', + 'outputdefinition', + 'results', + 'radaroverlay', + 'miscellaneous', 'private'] # }}} @@ -198,5 +197,5 @@ string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("levelset","[%s,%s]" % ("1x1",obj.levelset.__class__.__name__),"parameters for moving boundaries (level-set method)")) string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("calving","[%s,%s]" % ("1x1",obj.calving.__class__.__name__),"parameters for calving")) - string="%s\n%s" % (string,'%19s: %-22s -- %s' % ("love","[%s,%s]" % ("1x1",obj.love.__class__.__name__),"parameters for love solution")) + string="%s\n%s" % (string,'%19s: %-22s -- %s' % ("love","[%s,%s]" % ("1x1",obj.love.__class__.__name__),"parameters for love solution")) string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("autodiff","[%s,%s]" % ("1x1",obj.autodiff.__class__.__name__),"automatic differentiation parameters")) string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("inversion","[%s,%s]" % ("1x1",obj.inversion.__class__.__name__),"parameters for inverse methods")) @@ -293,5 +292,5 @@ field=getattr(md1,fieldi) fieldsize=np.shape(field) - if hasattr(field,'__dict__') and not m.ismember(fieldi,['results'])[0]: #recursive call + if hasattr(field,'__dict__') and not fieldi in ['results']: #recursive call object_fields=vars(field) for fieldj in object_fields: @@ -301,5 +300,5 @@ if len(fieldsize): #size = number of nodes * n - if fieldsize[0]==numberofvertices1: + if fieldsize[0]==numberofvertices1: setattr(getattr(md2,fieldi),fieldj,field[pos_node]) elif fieldsize[0]==numberofvertices1+1: @@ -311,5 +310,5 @@ if len(fieldsize): #size = number of nodes * n - if fieldsize[0]==numberofvertices1: + if fieldsize[0]==numberofvertices1: setattr(md2,fieldi,field[pos_node]) elif fieldsize[0]==numberofvertices1+1: @@ -362,11 +361,11 @@ #Edges - if m.strcmp(md.mesh.domaintype(),'2Dhorizontal'): + if md.mesh.domaintype()=='2Dhorizontal': if np.ndim(md2.mesh.edges)>1 and np.size(md2.mesh.edges,axis=1)>1: #do not use ~isnan because there are some np.nans... #renumber first two columns pos=np.nonzero(md2.mesh.edges[:,3]!=-1)[0] - md2.mesh.edges[: ,0]=Pnode[md2.mesh.edges[:,0]-1] - md2.mesh.edges[: ,1]=Pnode[md2.mesh.edges[:,1]-1] - md2.mesh.edges[: ,2]=Pelem[md2.mesh.edges[:,2]-1] + md2.mesh.edges[:,0]=Pnode[md2.mesh.edges[:,0]-1] + md2.mesh.edges[:,1]=Pnode[md2.mesh.edges[:,1]-1] + md2.mesh.edges[:,2]=Pelem[md2.mesh.edges[:,2]-1] md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3]-1] #remove edges when the 2 vertices are not in the domain. @@ -695,44 +694,69 @@ #dealing with the friction law #drag is limited to nodes that are on the bedrock. - if hasattr(md.friction,'coefficient'): md.friction.coefficient=project2d(md,md.friction.coefficient,1) + 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'): md.friction.p=project2d(md,md.friction.p,1) - if hasattr(md.friction,'q'): md.friction.q=project2d(md,md.friction.q,1) + if hasattr(md.friction,'p'): + md.friction.p=project2d(md,md.friction.p,1) + if hasattr(md.friction,'q'): + md.friction.q=project2d(md,md.friction.q,1) - if hasattr(md.friction,'coefficientcoulomb'): md.friction.coefficientcoulomb=project2d(md,md.friction.coefficientcoulomb,1) - if hasattr(md.friction,'C'): md.friction.C=project2d(md,md.friction.C,1) - if hasattr(md.friction,'As'): md.friction.As=project2d(md,md.friction.As,1) + if hasattr(md.friction,'coefficientcoulomb'): + md.friction.coefficientcoulomb=project2d(md,md.friction.coefficientcoulomb,1) + if hasattr(md.friction,'C'): + md.friction.C=project2d(md,md.friction.C,1) + 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'): md.friction.water_layer=project2d(md,md.friction.water_layer,1) - if hasattr(md.friction,'m'): md.friction.m=project2d(md,md.friction.m,1) + md.friction.effective_pressure=project2d(md,md.friction.effective_pressure,1) + if hasattr(md.friction,'water_layer'): + md.friction.water_layer=project2d(md,md.friction.water_layer,1) + 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 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 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 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 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 not np.isnan(md.smb.mass_balance).all(): md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers) #results - if not np.isnan(md.initialization.vx).all(): md.initialization.vx=DepthAverage(md,md.initialization.vx) - if not np.isnan(md.initialization.vy).all(): md.initialization.vy=DepthAverage(md,md.initialization.vy) - if not np.isnan(md.initialization.vz).all(): md.initialization.vz=DepthAverage(md,md.initialization.vz) - if not np.isnan(md.initialization.vel).all(): md.initialization.vel=DepthAverage(md,md.initialization.vel) - if not np.isnan(md.initialization.temperature).all(): md.initialization.temperature=DepthAverage(md,md.initialization.temperature) - 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(): md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1) - 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(): md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1) + if not np.isnan(md.initialization.vx).all(): + md.initialization.vx=DepthAverage(md,md.initialization.vx) + if not np.isnan(md.initialization.vy).all(): + md.initialization.vy=DepthAverage(md,md.initialization.vy) + if not np.isnan(md.initialization.vz).all(): + md.initialization.vz=DepthAverage(md,md.initialization.vz) + if not np.isnan(md.initialization.vel).all(): + md.initialization.vel=DepthAverage(md,md.initialization.vel) + if not np.isnan(md.initialization.temperature).all(): + md.initialization.temperature=DepthAverage(md,md.initialization.temperature) + 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(): + md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1) + 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(): + 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 @@ -760,6 +784,7 @@ md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers) md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers) - if not np.isnan(md.damage.spcdamage).all(): md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers-1) md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers-1) + if not np.isnan(md.damage.spcdamage).all(): + md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers-1) #materials @@ -787,6 +812,6 @@ md.geometry.thickness=project2d(md,md.geometry.thickness,1) md.geometry.base=project2d(md,md.geometry.base,1) - if isinstance(md.geometry.bed,np.ndarray): - md.geometry.bed=project2d(md,md.geometry.bed,1) + if isinstance(md.geometry.bed,np.ndarray): + md.geometry.bed=project2d(md,md.geometry.bed,1) md.mask.groundedice_levelset=project2d(md,md.mask.groundedice_levelset,1) md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1) @@ -815,10 +840,14 @@ mesh.numberofelements=md.mesh.numberofelements2d mesh.elements=md.mesh.elements2d - if not np.isnan(md.mesh.vertexonboundary).all(): mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1) - if not np.isnan(md.mesh.elementconnectivity).all(): mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1) + if not np.isnan(md.mesh.vertexonboundary).all(): + mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1) + 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: 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 md.mesh=mesh