Changeset 13525
- Timestamp:
- 10/03/12 15:56:53 (12 years ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/test/NightlyRun/test217.py
r13497 r13525 19 19 from parameterize import * 20 20 from setflowequation import * 21 from paterson import * 21 22 from solve import * 22 23 … … 27 28 md.cluster=generic('name',oshostname(),'np',3) 28 29 29 30 30 # redo the parameter file for this special shelf. 31 32 31 # constant thickness, constrained (vy=0) flow into an icefront, 33 34 32 # from 0 m/yr at the grounding line. 35 33 36 37 38 34 # tighten 39 40 md.diagnostic.restol=10^-4 41 35 md.diagnostic.restol=10**-4 42 36 43 37 # needed later 44 45 38 ymin=min(md.mesh.y) 46 39 ymax=max(md.mesh.y) … … 48 41 xmax=max(md.mesh.x) 49 42 50 51 43 di=md.materials.rho_ice/md.materials.rho_water 52 44 53 54 h=1000 55 md.geometry.thickness=h*ones(md.mesh.numberofvertices,1) 45 h=1000. 46 md.geometry.thickness=h*ones((md.mesh.numberofvertices,1)) 56 47 md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness 57 48 md.geometry.surface=md.geometry.bed+md.geometry.thickness 58 49 59 60 50 # Initial velocity and pressure 61 62 md.initialization.vx=zeros(md.mesh.numberofvertices,1) 63 md.initialization.vy=zeros(md.mesh.numberofvertices,1) 64 md.initialization.vz=zeros(md.mesh.numberofvertices,1) 65 md.initialization.pressure=zeros(md.mesh.numberofvertices,1) 66 51 md.initialization.vx=zeros((md.mesh.numberofvertices,1)) 52 md.initialization.vy=zeros((md.mesh.numberofvertices,1)) 53 md.initialization.vz=zeros((md.mesh.numberofvertices,1)) 54 md.initialization.pressure=zeros((md.mesh.numberofvertices,1)) 67 55 68 56 # Materials 69 70 md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1) 57 md.initialization.temperature=(273.-20.)*ones((md.mesh.numberofvertices,1)) 71 58 md.materials.rheology_B=paterson(md.initialization.temperature) 72 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1) 73 59 md.materials.rheology_n=3*ones((md.mesh.numberofelements,1)) 74 60 75 61 # Boundary conditions: 76 77 md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1) 78 md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1) 79 md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1) 80 62 md.diagnostic.spcvx=float(nan)*ones((md.mesh.numberofvertices,1)) 63 md.diagnostic.spcvy=float(nan)*ones((md.mesh.numberofvertices,1)) 64 md.diagnostic.spcvz=float(nan)*ones((md.mesh.numberofvertices,1)) 81 65 82 66 # constrain flanks to 0 normal velocity 83 84 pos=numpy.nonzero(md.mesh.x==xmin | md.mesh.x==xmax) 85 md.diagnostic.spcvx(pos)=0 86 md.diagnostic.spcvz(pos)=NaN 87 67 pos=numpy.nonzero(numpy.logical_or(md.mesh.x==xmin,md.mesh.x==xmax)) 68 md.diagnostic.spcvx[pos]=0 69 md.diagnostic.spcvz[pos]=float(nan) 88 70 89 71 # constrain grounding line to 0 velocity 90 91 72 pos=numpy.nonzero(md.mesh.y==ymin) 92 md.diagnostic.spcvx(pos)=0 93 md.diagnostic.spcvy(pos)=0 94 73 md.diagnostic.spcvx[pos]=0 74 md.diagnostic.spcvy[pos]=0 95 75 96 76 # icefront 97 98 nodeonicefront=zeros(md.mesh.numberofvertices,1) 77 nodeonicefront=zeros(md.mesh.numberofvertices) 99 78 pos=numpy.nonzero(md.mesh.y==ymax) 100 nodeonicefront(pos)=1 101 pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2))) diagnostic.icefront=md.mesh.segments(pos,:) 102 diagnostic.icefront=[diagnostic.icefront 1*md.mask.elementonfloatingice(diagnostic.icefront(:,end))] 79 nodeonicefront[pos]=1 80 pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0].astype(int)-1],nodeonicefront[md.mesh.segments[:,1].astype(int)-1]))[0] 81 diagnostic.icefront=md.mesh.segments[pos,:] 82 diagnostic.icefront=numpy.hstack((diagnostic.icefront,1.*md.mask.elementonfloatingice[diagnostic.icefront[:,-1].astype(int)-1].reshape(-1,1))) 103 83 md.diagnostic.icefront=diagnostic.icefront 104 105 84 106 85 md=solve(md,DiagnosticSolutionEnum()) 107 86 108 109 87 # create analytical solution: strain rate is constant = ((rho_ice*g*h)/4B)^3 (Paterson, 4th Edition, page 292. 110 111 88 # ey_c=(md.materials.rho_ice*md.constants.g*(1-di)*md.geometry.thickness./(4*md.materials.rheology_B)).^3; 112 113 89 # vy_c=ey_c.*md.mesh.y*md.constants.yts; 114 90 115 116 117 91 # Fields and tolerances to track changes 118 119 92 field_names =['Vy'] 120 93 field_tolerances=[1e-13]
Note:
See TracChangeset
for help on using the changeset viewer.