Changeset 13525

Timestamp:

10/03/12 15:56:53 (12 years ago)

Author:

jschierm

Message:

NEW: Working test217.py.

File:

• issm/trunk-jpl/test/NightlyRun/test217.py (moved) (moved from issm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test217.py ) (3 diffs)

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

@@ -19 +19 @@
 from parameterize import *
 from setflowequation import *
+from paterson import *
 from solve import *
 
@@ -27 +28 @@
 md.cluster=generic('name',oshostname(),'np',3)
 
-
 # redo the parameter file for this special shelf. 
-
 # constant thickness, constrained (vy=0) flow into an icefront, 
-
 # from 0 m/yr at the grounding line.
 
-
-
 # tighten
-
-md.diagnostic.restol=10^-4
-
+md.diagnostic.restol=10**-4
 
 # needed later
-
 ymin=min(md.mesh.y)
 ymax=max(md.mesh.y)
@@ -48 +41 @@
 xmax=max(md.mesh.x)
 
-
 di=md.materials.rho_ice/md.materials.rho_water
 
-
-h=1000
-md.geometry.thickness=h*ones(md.mesh.numberofvertices,1)
+h=1000.
+md.geometry.thickness=h*ones((md.mesh.numberofvertices,1))
 md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
 md.geometry.surface=md.geometry.bed+md.geometry.thickness
 
-
 # Initial velocity and pressure
-
-md.initialization.vx=zeros(md.mesh.numberofvertices,1)
-md.initialization.vy=zeros(md.mesh.numberofvertices,1)
-md.initialization.vz=zeros(md.mesh.numberofvertices,1)
-md.initialization.pressure=zeros(md.mesh.numberofvertices,1)
-
+md.initialization.vx=zeros((md.mesh.numberofvertices,1))
+md.initialization.vy=zeros((md.mesh.numberofvertices,1))
+md.initialization.vz=zeros((md.mesh.numberofvertices,1))
+md.initialization.pressure=zeros((md.mesh.numberofvertices,1))
 
 # Materials
-
-md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1)
+md.initialization.temperature=(273.-20.)*ones((md.mesh.numberofvertices,1))
 md.materials.rheology_B=paterson(md.initialization.temperature)
-md.materials.rheology_n=3*ones(md.mesh.numberofelements,1)
-
+md.materials.rheology_n=3*ones((md.mesh.numberofelements,1))
 
 # Boundary conditions:
-
-md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1)
-md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1)
-md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1)
-
+md.diagnostic.spcvx=float(nan)*ones((md.mesh.numberofvertices,1))
+md.diagnostic.spcvy=float(nan)*ones((md.mesh.numberofvertices,1))
+md.diagnostic.spcvz=float(nan)*ones((md.mesh.numberofvertices,1))
 
 # constrain flanks to 0 normal velocity
-
-pos=numpy.nonzero(md.mesh.x==xmin | md.mesh.x==xmax)
-md.diagnostic.spcvx(pos)=0
-md.diagnostic.spcvz(pos)=NaN
-
+pos=numpy.nonzero(numpy.logical_or(md.mesh.x==xmin,md.mesh.x==xmax))
+md.diagnostic.spcvx[pos]=0
+md.diagnostic.spcvz[pos]=float(nan)
 
 # constrain grounding line to 0 velocity
-
 pos=numpy.nonzero(md.mesh.y==ymin)
-md.diagnostic.spcvx(pos)=0
-md.diagnostic.spcvy(pos)=0
-
+md.diagnostic.spcvx[pos]=0
+md.diagnostic.spcvy[pos]=0
 
 # icefront
-
-nodeonicefront=zeros(md.mesh.numberofvertices,1)
+nodeonicefront=zeros(md.mesh.numberofvertices)
 pos=numpy.nonzero(md.mesh.y==ymax)
-nodeonicefront(pos)=1
-pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2))) diagnostic.icefront=md.mesh.segments(pos,:)
-diagnostic.icefront=[diagnostic.icefront 1*md.mask.elementonfloatingice(diagnostic.icefront(:,end))]
+nodeonicefront[pos]=1
+pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0].astype(int)-1],nodeonicefront[md.mesh.segments[:,1].astype(int)-1]))[0]
+diagnostic.icefront=md.mesh.segments[pos,:]
+diagnostic.icefront=numpy.hstack((diagnostic.icefront,1.*md.mask.elementonfloatingice[diagnostic.icefront[:,-1].astype(int)-1].reshape(-1,1)))
 md.diagnostic.icefront=diagnostic.icefront
-
 
 md=solve(md,DiagnosticSolutionEnum())
 
-
 # create analytical solution: strain rate is constant = ((rho_ice*g*h)/4B)^3 (Paterson, 4th Edition, page 292.
-
 # ey_c=(md.materials.rho_ice*md.constants.g*(1-di)*md.geometry.thickness./(4*md.materials.rheology_B)).^3;
-
 # vy_c=ey_c.*md.mesh.y*md.constants.yts;
 
-
-
 # Fields and tolerances to track changes
-
 field_names     =['Vy']
 field_tolerances=[1e-13]