Changeset 13525


Ignore:
Timestamp:
10/03/12 15:56:53 (12 years ago)
Author:
jschierm
Message:

NEW: Working test217.py.

File:
1 moved

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/test/NightlyRun/test217.py

    r13497 r13525  
    1919from parameterize import *
    2020from setflowequation import *
     21from paterson import *
    2122from solve import *
    2223
     
    2728md.cluster=generic('name',oshostname(),'np',3)
    2829
    29 
    3030# redo the parameter file for this special shelf.
    31 
    3231# constant thickness, constrained (vy=0) flow into an icefront,
    33 
    3432# from 0 m/yr at the grounding line.
    3533
    36 
    37 
    3834# tighten
    39 
    40 md.diagnostic.restol=10^-4
    41 
     35md.diagnostic.restol=10**-4
    4236
    4337# needed later
    44 
    4538ymin=min(md.mesh.y)
    4639ymax=max(md.mesh.y)
     
    4841xmax=max(md.mesh.x)
    4942
    50 
    5143di=md.materials.rho_ice/md.materials.rho_water
    5244
    53 
    54 h=1000
    55 md.geometry.thickness=h*ones(md.mesh.numberofvertices,1)
     45h=1000.
     46md.geometry.thickness=h*ones((md.mesh.numberofvertices,1))
    5647md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
    5748md.geometry.surface=md.geometry.bed+md.geometry.thickness
    5849
    59 
    6050# 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 
     51md.initialization.vx=zeros((md.mesh.numberofvertices,1))
     52md.initialization.vy=zeros((md.mesh.numberofvertices,1))
     53md.initialization.vz=zeros((md.mesh.numberofvertices,1))
     54md.initialization.pressure=zeros((md.mesh.numberofvertices,1))
    6755
    6856# Materials
    69 
    70 md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1)
     57md.initialization.temperature=(273.-20.)*ones((md.mesh.numberofvertices,1))
    7158md.materials.rheology_B=paterson(md.initialization.temperature)
    72 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1)
    73 
     59md.materials.rheology_n=3*ones((md.mesh.numberofelements,1))
    7460
    7561# 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 
     62md.diagnostic.spcvx=float(nan)*ones((md.mesh.numberofvertices,1))
     63md.diagnostic.spcvy=float(nan)*ones((md.mesh.numberofvertices,1))
     64md.diagnostic.spcvz=float(nan)*ones((md.mesh.numberofvertices,1))
    8165
    8266# 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 
     67pos=numpy.nonzero(numpy.logical_or(md.mesh.x==xmin,md.mesh.x==xmax))
     68md.diagnostic.spcvx[pos]=0
     69md.diagnostic.spcvz[pos]=float(nan)
    8870
    8971# constrain grounding line to 0 velocity
    90 
    9172pos=numpy.nonzero(md.mesh.y==ymin)
    92 md.diagnostic.spcvx(pos)=0
    93 md.diagnostic.spcvy(pos)=0
    94 
     73md.diagnostic.spcvx[pos]=0
     74md.diagnostic.spcvy[pos]=0
    9575
    9676# icefront
    97 
    98 nodeonicefront=zeros(md.mesh.numberofvertices,1)
     77nodeonicefront=zeros(md.mesh.numberofvertices)
    9978pos=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))]
     79nodeonicefront[pos]=1
     80pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0].astype(int)-1],nodeonicefront[md.mesh.segments[:,1].astype(int)-1]))[0]
     81diagnostic.icefront=md.mesh.segments[pos,:]
     82diagnostic.icefront=numpy.hstack((diagnostic.icefront,1.*md.mask.elementonfloatingice[diagnostic.icefront[:,-1].astype(int)-1].reshape(-1,1)))
    10383md.diagnostic.icefront=diagnostic.icefront
    104 
    10584
    10685md=solve(md,DiagnosticSolutionEnum())
    10786
    108 
    10987# create analytical solution: strain rate is constant = ((rho_ice*g*h)/4B)^3 (Paterson, 4th Edition, page 292.
    110 
    11188# ey_c=(md.materials.rho_ice*md.constants.g*(1-di)*md.geometry.thickness./(4*md.materials.rheology_B)).^3;
    112 
    11389# vy_c=ey_c.*md.mesh.y*md.constants.yts;
    11490
    115 
    116 
    11791# Fields and tolerances to track changes
    118 
    11992field_names     =['Vy']
    12093field_tolerances=[1e-13]
Note: See TracChangeset for help on using the changeset viewer.