Python PISM类代码示例

    def solve(self):
        """Solve the SSA by calling the underlying PISM :cpp:class:`SSA`'s
        :cpp:member:`update` method. Returns the solution vector (owned by
        self.ssa, but you should not need to know about ownership).

        """
        vecs = self.modeldata.vecs

        # make sure vecs is locked!
        self.ssa.init()
        if vecs.has('vel_bc'):
            self.ssa.set_boundary_conditions(vecs.bc_mask, vecs.vel_bc)

        melange_back_pressure = PISM.IceModelVec2S()
        melange_back_pressure.create(self.grid, "melange_back_pressure", PISM.WITHOUT_GHOSTS)
        melange_back_pressure.set_attrs("diagnostic",
                                        "melange back pressure fraction", "1", "")

        PISM.verbPrintf(2, self.grid.com, "* Solving the SSA stress balance ...\n")

        fast = False

        self.ssa.update(fast, melange_back_pressure)

        return self.ssa.velocity()

    def write(self, filename):
        """Saves all of :attr:`modeldata`'s vecs (and the solution) to an
        output file."""
        grid = self.grid
        vecs = self.modeldata.vecs

        pio = PISM.PIO(grid.com, "netcdf3")
        pio.open(filename, PISM.PISM_READWRITE_MOVE)
        PISM.define_time(pio, grid.ctx().config().get_string("time_dimension_name"),
                         grid.ctx().config().get_string("calendar"),
                         grid.ctx().time().units_string(),
                         grid.ctx().unit_system())
        PISM.append_time(pio, grid.ctx().config().get_string("time_dimension_name"), 0.0)
        pio.close()

        # Save time & command line
        PISM.util.writeProvenance(filename)

        vecs.writeall(filename)

        vel_ssa = self.ssa.velocity()
        vel_ssa.write(filename)

        sys = self.grid.ctx().unit_system()

        velbar_mag = model.createCBarVec(self.grid)
        velbar_mag.set_to_magnitude(vel_ssa)
        velbar_mag.mask_by(vecs.thk, PISM.convert(sys, -0.01, "m/year", "m/second"))
        velbar_mag.write(filename)

    def solve(self):
        """Solve the SSA by calling the underlying PISM :cpp:class:`SSA`'s
        :cpp:member:`update` method. Returns the solution vector (owned by
        self.ssa, but you should not need to know about ownership).

        """
        vecs = self.modeldata.vecs

        # make sure vecs is locked!
        self.ssa.init()

        melange_back_pressure = PISM.IceModelVec2S()
        melange_back_pressure.create(self.grid, "melange_back_pressure", PISM.WITHOUT_GHOSTS)
        melange_back_pressure.set_attrs("diagnostic",
                                        "melange back pressure fraction", "1", "")
        melange_back_pressure.set(0.0)

        PISM.verbPrintf(2, self.grid.com, "* Solving the SSA stress balance ...\n")

        full_update = True

        inputs                       = PISM.StressBalanceInputs()
        inputs.melange_back_pressure = melange_back_pressure
        inputs.geometry              = self.geometry
        inputs.enthalpy              = vecs.enthalpy
        inputs.basal_yield_stress    = vecs.tauc
        if vecs.has('vel_bc'):
            inputs.bc_mask   = vecs.bc_mask
            inputs.bc_values = vecs.vel_bc

        self.ssa.update(inputs, full_update)

        return self.ssa.velocity()

def grounded_cell_fraction_test():

    # allocation
    grid = allocate_grid(ctx)

    ice_thickness, bed_topography, surface, mask, gl_mask, gl_mask_x, gl_mask_y, _ = allocate_storage(grid)

    bed_topography.set(0.0)

    # initialization
    sea_level = 500.0
    for L in [0.0, 0.25, 0.5, 0.75, 1.0]:
        init(mu, L, sea_level, ice_thickness, "box")

        compute_mask(sea_level, bed_topography, ice_thickness, mask, surface)

        # computation of gl_mask
        PISM.compute_grounded_cell_fraction(ice_density, ocean_density, sea_level,
                                            ice_thickness, bed_topography, mask, gl_mask,
                                            gl_mask_x, gl_mask_y)

        # inspection / comparison
        print("L = %f" % L)
        print_vec(mask)
        print_vec(gl_mask_x)
        print_vec(gl_mask_y)
        print_vec(gl_mask)

  def report_riggs(self):
    grid = self.grid
    r=PISM.VerbPrintf(grid.com,verbosity=2)
    r.println("comparing to RIGGS data in %s ...\n",self.riggs_file); 

    riggsvars = [ "riggslat", "riggslon", "riggsmag", "riggsu", "riggsv"]
    (latdata,longdata,magdata,udata,vdata) = \
         (PISM.Timeseries(grid,varname,"count") for varname in riggsvars)
    riggsvars = [latdata, longdata, magdata, udata, vdata]

    for v in [magdata, udata, vdata]:
      v.set_units("m year-1", "")

    for v in riggsvars:
      v.read(self.riggs_file)

    length = latdata.length();

    vel_ssa = self.solver.solution();
    clat = self.latitude; clon = self.longitude; mask = self.solver.ice_mask;

    secpera=PISM.secpera
    with PISM.util.Access([clat,clon,mask,vel_ssa]):
      goodptcount = 0.0; ChiSqr = 0.0;    
      for k in xrange(length):
        (lat,lon,mag,u,v) = [v[k] for v in riggsvars]
        r.printlnv(4," RIGGS[%3d]: lat = %7.3f, lon = %7.3f, mag = %7.2f, u = %7.2f, v = %7.2f",
                          k,lat,lon,mag,u,v)
        origdlat = (-5.42445 - (-12.3325)) / 110.0;
        lowlat = -12.3325 - origdlat * 46.0;
        dlat = (-5.42445 - lowlat) / (float) (grid.My - 1);        
        lowlon = -5.26168;
        dlon = (3.72207 - lowlon) / (float) (grid.Mx - 1);
        cj = int( math.floor((lat - lowlat) / dlat) )
        ci = int( math.floor((lon - lowlon) / dlon) )
      
        if ((ci >= grid.xs) and (ci < grid.xs+grid.xm) and (cj >= grid.ys) and (cj < grid.ys+grid.ym)):
          vel = vel_ssa[ci,cj]
          cu = secpera * vel.u
          cv = secpera * vel.v
          cmag = math.sqrt(cu*cu + cv*cv)
          PISM.verbPrintf(4,PETSc.COMM_SELF,
                          " PISM%d[%3d]: lat = %7.3f, lon = %7.3f, mag = %7.2f, u = %7.2f, v = %7.2f\n",
                          grid.rank,k,clat[ci,cj],clon[ci,cj],cmag,cu,cv)
          if mask[ci,cj] == PISM.MASK_FLOATING:
            goodptcount += 1.0;
            ChiSqr += (u-cu)*(u-cu)+(v-cv)*(v-cv)
    # end with

    ChiSqr = ChiSqr / (30.0*30.0) # see page 48 of MacAyeal et al
    g_goodptcount = PISM.globalSum(goodptcount,grid.com)
    g_ChiSqr      = PISM.globalSum(ChiSqr, grid.com)
    r.printlnv( 4, """number of RIGGS data points = %d
number of RIGGS points in computed ice shelf region = %8.2f""", length , g_goodptcount);
    r.println("Chi^2 statistic for computed results compared to RIGGS is %10.3f",
                      g_ChiSqr * (156.0 / g_goodptcount))

def pause(message_in=None, message_out=None):
    """Prints a message and waits for a key press.

    :param message_in: Message to display before waiting.
    :param message_out: Message to display after waiting."""
    com = PISM.Context().com
    if not message_in is None:
        PISM.verbPrintf(1, com, message_in + "\n")
    _ = getkey()
    if not message_out is None:
        PISM.verbPrintf(1, com, message_out + "\n")

def pism_pause(message_in=None, message_out=None):
    """Implements a ``siple`` pause callback appropriate for parallel runs."""
    import sys
    import os
    fd = sys.stdin.fileno()
    com = PISM.Context().com
    if os.isatty(fd):
        return siple.reporting.std_pause(message_in, message_out)
    if not message_in is None:
        PISM.verbPrintf(1, com, message_in + "\n")
    _ = sys.stdin.read(1)
    if not message_out is None:
        PISM.verbPrintf(1, com, message_out + "\n")

def writeProvenance(outfile, message=None):
    """Saves the time and command line arguments (or the provided `message`) to
    the ``history`` attribute of the :file:`.nc` file `outfile`"""

    com = PISM.Context().com

    ds = PISM.PIO(com, "netcdf3", outfile, PISM.PISM_READWRITE)

    if message is None:
        message = PISM.timestamp(com) + ": " + PISM.args_string()
    ds.append_history(message)
    ds.put_att_text("PISM_GLOBAL", "source", "PISM " + PISM.PISM_Revision)

    ds.close()

def prepare_output(filename, append_time=True):
    "Create an output file and prepare it for writing."
    ctx = PISM.Context()
    output = PISM.PIO(ctx.com, ctx.config.get_string("output.format"),
                      filename, PISM.PISM_READWRITE_MOVE)
    PISM.define_time(output,
                     ctx.config.get_string("time.dimension_name"),
                     ctx.config.get_string("time.calendar"),
                     ctx.time.units_string(),
                     ctx.unit_system)
    if append_time:
        PISM.append_time(output,
                         ctx.config.get_string("time.dimension_name"),
                         ctx.time.current())

    return output

  def updateYieldStress(self,mask,thickness,Hmelt,bmr,tillphi,tauc):
    config = PISM.global_config()
    till_pw_fraction = self.till_pw_fraction#config.get("till_pw_fraction")
    till_c_0 = self.till_c_0#config.get("till_c_0") * 1e3 # convert from kPa to Pa
    hmelt_max = self.hmelt_max#config.get("hmelt_max");

    rho_g = self.rho_g


    with PISM.util.Access(nocomm=[mask,tauc,thickness,Hmelt,bmr,tillphi]):
      mq = PISM.MaskQuery(mask)
      GHOSTS = self.grid.max_stencil_width;
      for (i,j) in self.grid.points_with_ghosts(nGhosts = GHOSTS):
        if mq.floating_ice(i,j):
          tauc[i,j] = 0
          continue

        H_ij = thickness[i,j]
        if H_ij == 0:
          tauc[i,j] = 1000.0e3;  #// large yield stress of 1000 kPa = 10 bar if no ice
        else: # grounded and there is some ice
          p_over = rho_g * H_ij
          p_w    = self.getBasalWaterPressure( H_ij,
                                         Hmelt[i,j],bmr[i,j],till_pw_fraction, 
                                         hmelt_max)
          N = p_over - p_w #  effective pressure on till
          tauc[i,j] = till_c_0 + N * math.tan((math.pi/180.0) * tillphi[i,j])

def create2dVelocityVec(grid, name="", desc="", intent="", ghost_type=PISM.WITH_GHOSTS, stencil_width=None):
    """Returns an :cpp:class:`IceModelVec2V` with attributes setup for horizontal velocity data.

    :param grid: The grid to associate with the vector.
    :param name: The intrinsic name to give the vector.

    :param ghost_type: One of ``PISM.WITH_GHOSTS`` or
                       ``PISM.WITHOUT_GHOSTS`` indicating if the vector
                       is ghosted.

    :param stencil_width: The size of the ghost stencil. Ignored if
                          ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``. If
                          ``stencil_width`` is ``None`` and the vector
                          is ghosted, the grid's maximum stencil width
                          is used.

    """
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    vel = PISM.IceModelVec2V()
    vel.create(grid, name, ghost_type, stencil_width)
    vel.set_attrs(intent, "%s%s" % ("X-component of the ", desc), "m s-1", "", 0)
    vel.set_attrs(intent, "%s%s" % ("Y-component of the ", desc), "m s-1", "", 1)
    vel.metadata(0).set_string("glaciological_units", "m year-1")
    vel.metadata(1).set_string("glaciological_units", "m year-1")
    vel.write_in_glaciological_units = True
    sys = grid.ctx().unit_system()
    huge_vel = PISM.convert(sys, 1e10, "m/year", "m/second")
    attrs = [("valid_min", -huge_vel), ("valid_max", huge_vel), ("_FillValue", 2 * huge_vel)]
    for a in attrs:
        for component in [0, 1]:
            vel.metadata(component).set_double(a[0], a[1])
    vel.set(2 * huge_vel)
    return vel

def flowlaw_test():
    ctx = PISM.context_from_options(PISM.PETSc.COMM_WORLD, "flowlaw_test")
    EC = ctx.enthalpy_converter()
    ff = PISM.FlowLawFactory("sia_", ctx.config(), EC)
    law = ff.create()

    TpaC = [-30, -5, 0, 0]
    depth = 2000
    gs = 1e-3
    omega = [0.0, 0.0, 0.0, 0.005]
    sigma = [1e4, 5e4, 1e5, 1.5e5]

    p = EC.pressure(depth)
    Tm = EC.melting_temperature(p)

    print "flow law:   \"%s\"" % law.name()
    print "pressure = %9.3e Pa = (hydrostatic at depth %7.2f m)" % (p, depth)
    print "flowtable:"
    print "  (dev stress)   (abs temp) (liq frac) =   (flow)"

    for i in range(4):
        for j in range(4):
            T = Tm + TpaC[j]
            E = EC.enthalpy(T, omega[j], p)
            flowcoeff = law.flow(sigma[i], E, p, gs)
            print "    %10.2e   %10.3f  %9.3f = %10.6e" % (sigma[i], T, omega[j], flowcoeff)

def add_disc_load(ice_thickness, radius, thickness):
    "Add a disc load with a given radius and thickness."
    grid = ice_thickness.grid()

    with PISM.vec.Access(nocomm=ice_thickness):
        for (i, j) in grid.points():
            r = PISM.radius(grid, i, j)
            if r <= disc_radius:
                ice_thickness[i, j] = disc_thickness

 def __call__(self, message, verbosity):
     """Saves the message to our internal log string and writes the string out to the file."""
     if self.rank == 0 and PISM.getVerbosityLevel() >= verbosity:
         timestamp = time.strftime('%Y-%m-%d %H:%M:%S')
         self.log = "%s%s: %s" % (self.log, timestamp, message)
         d = PISM.netCDF.Dataset(self.filename, 'a')
         d.__setattr__(self.attr, self.log)
         d.close()
     self.com.barrier()

  def initPhysics(self):
    config = self.config
    self.basal = PISM.IceBasalResistancePlasticLaw(
           config.get("plastic_regularization") / PISM.secpera,
           config.get_flag("do_pseudo_plastic_till"),
           config.get("pseudo_plastic_q"),
           config.get("pseudo_plastic_uthreshold") / PISM.secpera);

    if PISM.optionsIsSet("-ssa_glen"):
      self.ice = PISM.CustomGlenIce(self.grid.com,"",config)
      B_schoof = 3.7e8;     # Pa s^{1/3}; hardness 
      self.ice.setHardness(B_schoof)
    else:
      self.ice =  PISM.GPBLDIce(self.grid.com, "", config)
    self.ice.setFromOptions()
    
    self.enthalpyconverter = PISM.EnthalpyConverter(config)
    if PISM.getVerbosityLevel() >3:
      self.enthalpyconverter.viewConstants(PETSc.Viewer.STDOUT())