Python PISM.convert方法代码示例

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import convert [as 别名]
    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)

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import convert [as 别名]
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

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import convert [as 别名]
    def _initSSA(self):
        # Test J has a viscosity that is independent of velocity.  So we force a
        # constant viscosity by settting the strength_extension
        # thickness larger than the given ice thickness. (max = 770m).

        sys = self.grid.ctx().unit_system()
        nu0 = PISM.convert(sys, 30.0, "MPa year", "Pa s")
        H0 = 500.0                        # 500 m typical thickness

        ssa = self.ssa
        ssa.strength_extension.set_notional_strength(nu0 * H0)
        ssa.strength_extension.set_min_thickness(800.)

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import convert [as 别名]
def ssa_trivial_test():
    "Test the SSA solver using a trivial setup."

    context = PISM.Context()
    unit_system = context.unit_system

    L = 50.e3  # // 50km half-width
    H0 = 500  # // m
    dhdx = 0.005  # // pure number, slope of surface & bed
    nu0 = PISM.convert(unit_system, 30.0, "MPa year", "Pa s")
    tauc0 = 1.e4  # // 1kPa

    class TrivialSSARun(PISM.ssa.SSAExactTestCase):
        def _initGrid(self):
            self.grid = PISM.IceGrid.Shallow(PISM.Context().ctx, L, L, 0, 0,
                                             self.Mx, self.My, PISM.NONE)

        def _initPhysics(self):
            self.modeldata.setPhysics(context.enthalpy_converter)

        def _initSSACoefficients(self):
            self._allocStdSSACoefficients()
            self._allocateBCs()

            vecs = self.modeldata.vecs

            vecs.land_ice_thickness.set(H0)
            vecs.surface_altitude.set(H0)
            vecs.bedrock_altitude.set(0.0)
            vecs.tauc.set(tauc0)

            # zero Dirichler B.C. everywhere
            vecs.vel_bc.set(0.0)
            vecs.bc_mask.set(1.0)

        def _initSSA(self):
            # The following ensure that the strength extension is used everywhere
            se = self.ssa.strength_extension
            se.set_notional_strength(nu0 * H0)
            se.set_min_thickness(4000 * 10)

            # For the benefit of SSAFD on a non-periodic grid
            self.config.set_boolean("compute_surf_grad_inward_ssa", True)

        def exactSolution(self, i, j, x, y):
            return [0, 0]

    Mx = 11
    My = 11
    test_case = TrivialSSARun(Mx, My)
    test_case.run("ssa_trivial.nc")

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import convert [as 别名]
def pism_context_test():
    "Test creating and using a C++-level Context"

    com = PISM.PETSc.COMM_WORLD
    system = PISM.UnitSystem("")

    logger = PISM.Logger(com, 2)

    config = PISM.DefaultConfig(com, "pism_config", "-config", system)
    config.init_with_default(logger)

    EC = PISM.EnthalpyConverter(config)

    time = PISM.Time(config, "360_day", system)

    ctx = PISM.cpp.Context(com, system, config, EC, time, logger, "greenland")

    print ctx.com().Get_size()
    print ctx.config().get_double("standard_gravity")
    print ctx.enthalpy_converter().L(273.15)
    print ctx.time().current()
    print PISM.convert(ctx.unit_system(), 1, "km", "m")
    print ctx.prefix()

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import convert [as 别名]
def po_constant_test():
    """Test that the basal melt rate computed by ocean::Constant is the
    same regardless of whether it is set using
    ocean_sub_shelf_heat_flux_into_ice or the command-line option."""

    grid = create_dummy_grid()
    config = grid.ctx().config()

    L = config.get_double("water_latent_heat_fusion")
    rho = config.get_double("ice_density")

    # prescribe a heat flux that corresponds to a mass flux which is
    # an integer multiple of m / year so that we can easily specify it
    # using a command-line option
    M = PISM.convert(grid.ctx().unit_system(), 1, "m / year", "m / second")
    Q_default = config.get_double("ocean_sub_shelf_heat_flux_into_ice")
    Q = M * L * rho
    config.set_double("ocean_sub_shelf_heat_flux_into_ice", Q)

    # without the command-line option
    ocean_constant = PISM.OceanConstant(grid)
    ocean_constant.init()
    mass_flux_1 = PISM.model.createIceThicknessVec(grid)
    ocean_constant.shelf_base_mass_flux(mass_flux_1)

    # reset Q
    config.set_double("ocean_sub_shelf_heat_flux_into_ice", Q_default)

    # with the command-line option
    o = PISM.PETSc.Options()
    o.setValue("-shelf_base_melt_rate", 1.0)

    ocean_constant = PISM.OceanConstant(grid)
    ocean_constant.init()
    mass_flux_2 = PISM.model.createIceThicknessVec(grid)
    ocean_constant.shelf_base_mass_flux(mass_flux_2)

    import numpy as np
    with PISM.vec.Access(nocomm=[mass_flux_1, mass_flux_2]):
        assert np.fabs(mass_flux_1[0, 0] - M * rho) < 1e-16
        assert np.fabs(mass_flux_2[0, 0] - M * rho) < 1e-16

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import convert [as 别名]
    def rangeVector(self):
        """Constructs a brand new vector from the range (i.e. state) vector space"""
        v = PISM.IceModelVec2V()
        v.create(self.grid, "", True, WIDE_STENCIL)

        # Add appropriate meta data.
        intent = "?inverse?"  # FIXME
        desc = "SSA velocity computed by inversion"
        v.set_attrs(intent, "%s%s" % ("X-component of the ", desc), "m s-1", "", 0)
        v.set_attrs(intent, "%s%s" % ("Y-component of the ", desc), "m s-1", "", 1)
        v.metadata(0).set_string("glaciological_units", "m year-1")
        v.metadata(1).set_string("glaciological_units", "m year-1")
        v.write_in_glaciological_units = True
        sys = self.grid.ctx().unit_system()
        huge_vel = PISM.convert(sys, 1e6, "m/year", "m/second")
        attrs = [("valid_min", -huge_vel), ("valid_max", huge_vel), ("_FillValue", 2 * huge_vel)]
        for a in attrs:
            for component in range(2):
                v.metadata(component).set_double(a[0], a[1])

        return PISMLocalVector(v)

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import convert [as 别名]
    def __call__(self, inverse_solver, count, data):

        if not self.l2_weight_init:
            vecs = inverse_solver.ssarun.modeldata.vecs
            if vecs.has('vel_misfit_weight'):
                self.l2_weight = self.toproczero(vecs.vel_misfit_weight)
            self.l2_weight_init = True

        method = inverse_solver.method

        r = self.toproczero(data.residual)
        Td = None
        if data.has_key('T_zeta_step'):
            Td = self.toproczero(data.T_zeta_step)
        TStarR = None
        if data.has_key('TStar_residual'):
            TStarR = self.toproczero(data.TStar_residual)
        d = None
        if data.has_key('zeta_step'):
            d = self.toproczero(data.zeta_step)
        zeta = self.toproczero(data.zeta)

        sys = self.grid.ctx().unit_system()
        secpera = PISM.convert(sys, 1.0, "year", "second")

        if self.grid.rank() == 0:
            import matplotlib.pyplot as pp

            pp.figure(self.figure())

            l2_weight = self.l2_weight

            pp.clf()

            V = self.Vmax

            pp.subplot(2, 3, 1)
            if l2_weight is not None:
                rx = l2_weight * r[0, :, :] * secpera
            else:
                rx = r[0, :, :] * secpera
            rx = np.maximum(rx, -V)
            rx = np.minimum(rx, V)
            pp.imshow(rx, origin='lower', interpolation='nearest')
            pp.colorbar()
            pp.title('r_x')
            pp.jet()

            pp.subplot(2, 3, 4)
            if l2_weight is not None:
                ry = l2_weight * r[1, :, :] * secpera
            else:
                ry = r[1, :, :] * secpera
            ry = np.maximum(ry, -V)
            ry = np.minimum(ry, V)
            pp.imshow(ry, origin='lower', interpolation='nearest')
            pp.colorbar()
            pp.title('r_y')
            pp.jet()

            if method == 'ign':
                pp.subplot(2, 3, 2)
                Tdx = Td[0, :, :] * secpera
                pp.imshow(Tdx, origin='lower', interpolation='nearest')
                pp.colorbar()
                pp.title('Td_x')
                pp.jet()

                pp.subplot(2, 3, 5)
                Tdy = Td[1, :, :] * secpera
                pp.imshow(Tdy, origin='lower', interpolation='nearest')
                pp.colorbar()
                pp.title('Td_y')
                pp.jet()
            elif method == 'sd' or method == 'nlcg':
                pp.subplot(2, 3, 2)
                pp.imshow(TStarR, origin='lower', interpolation='nearest')
                pp.colorbar()
                pp.title('TStarR')
                pp.jet()

            if d is not None:
                d *= -1
                pp.subplot(2, 3, 3)
                pp.imshow(d, origin='lower', interpolation='nearest')

                # colorbar does a divide by zero if 'd' is all zero,
                # as it will be at the start of iteration zero.
                # The warning message is a distraction, so we suppress it.
                import warnings
                with warnings.catch_warnings():
                    warnings.simplefilter("ignore")
                    pp.colorbar()
                    pp.jet()

                pp.title('-zeta_step')

            pp.subplot(2, 3, 6)
            pp.imshow(zeta, origin='lower', interpolation='nearest')
            pp.colorbar()
#.........这里部分代码省略.........

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import convert [as 别名]

#.........这里部分代码省略.........
        design = createDesignVec(grid, design_var)
        design.regrid(input_filename, True)
        design_prior.copy_from(design)
        vecs.add(design_prior, writing=True)

    if using_zeta_fixed_mask:
        if PISM.util.fileHasVariable(inv_data_filename, "zeta_fixed_mask"):
            zeta_fixed_mask = PISM.model.createZetaFixedMaskVec(grid)
            zeta_fixed_mask.regrid(inv_data_filename)
            vecs.add(zeta_fixed_mask)
        else:
            if design_var == 'tauc':
                logMessage("  Computing 'zeta_fixed_mask' (i.e. locations where design variable '%s' has a fixed value).\n" % design_var)
                zeta_fixed_mask = PISM.model.createZetaFixedMaskVec(grid)
                zeta_fixed_mask.set(1)
                mask = vecs.mask
                with PISM.vec.Access(comm=zeta_fixed_mask, nocomm=mask):
                    mq = PISM.MaskQuery(mask)
                    for (i, j) in grid.points():
                        if mq.grounded_ice(i, j):
                            zeta_fixed_mask[i, j] = 0
                vecs.add(zeta_fixed_mask)

                adjustTauc(vecs.mask, design_prior)
            elif design_var == 'hardav':
                PISM.logging.logPrattle("Skipping 'zeta_fixed_mask' for design variable 'hardav'; no natural locations to fix its value.")
                pass
            else:
                raise NotImplementedError("Unable to build 'zeta_fixed_mask' for design variable %s.", design_var)

    # Convert design_prior -> zeta_prior
    zeta_prior = PISM.IceModelVec2S()
    zeta_prior.create(grid, "zeta_prior", PISM.WITH_GHOSTS, WIDE_STENCIL)
    design_param.convertFromDesignVariable(design_prior, zeta_prior)
    vecs.add(zeta_prior, writing=True)

    # Determine the initial guess for zeta.  If we are restarting, load it from
    # the output file.  Otherwise, if 'zeta_inv' is in the inverse data file, use it.
    # If none of the above, copy from 'zeta_prior'.
    zeta = PISM.IceModelVec2S()
    zeta.create(grid, "zeta_inv", PISM.WITH_GHOSTS, WIDE_STENCIL)
    zeta.set_attrs("diagnostic", "zeta_inv", "1", "zeta_inv")
    if do_restart:
        # Just to be sure, verify that we have a 'zeta_inv' in the output file.
        if not PISM.util.fileHasVariable(output_filename, 'zeta_inv'):
            PISM.verbPrintf(1, com, "Unable to restart computation: file %s is missing variable 'zeta_inv'", output_filename)
            exit(1)
        PISM.logging.logMessage("  Inversion starting from 'zeta_inv' found in %s\n" % output_filename)
        zeta.regrid(output_filename, True)

    elif PISM.util.fileHasVariable(inv_data_filename, 'zeta_inv'):
        PISM.logging.logMessage("  Inversion starting from 'zeta_inv' found in %s\n" % inv_data_filename)
        zeta.regrid(inv_data_filename, True)
    else:
        zeta.copy_from(zeta_prior)

    vel_ssa_observed = None
    vel_ssa_observed = PISM.model.create2dVelocityVec(grid, '_ssa_observed', stencil_width=2)
    if PISM.util.fileHasVariable(inv_data_filename, "u_ssa_observed"):
        vel_ssa_observed.regrid(inv_data_filename, True)
        vecs.add(vel_ssa_observed, writing=saving_inv_data)
    else:
        if not PISM.util.fileHasVariable(inv_data_filename, "u_surface_observed"):
            PISM.verbPrintf(1, context.com, "Neither u/v_ssa_observed nor u/v_surface_observed is available in %s.\nAt least one must be specified.\n" % inv_data_filename)
            exit(1)
        vel_surface_observed = PISM.model.create2dVelocityVec(grid, '_surface_observed', stencil_width=2)

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import convert [as 别名]
  Does a time-independent calculation.  Does not run IceModel or a derived
  class thereof. Uses the van der Veen flow-line shelf geometry. Also may be
  used in a PISM software (regression) test.
"""

usage = \
    """
usage of SSA_TEST_CFBC:
  run ssa_test_cfbc -Mx <number> -My <number>
"""

context = PISM.Context()
unit_system = context.unit_system

H0 = 600.          # meters
V0 = PISM.convert(unit_system, 300, "m/year", "m/second")
C = 2.45e-18     # "typical constant ice parameter"
T = 400          # time used to compute the calving front location

Q0 = V0 * H0
Hc1 = 4. * C / Q0
Hc2 = 1. / (H0 ** 4)


def H_exact(x):
    return (Hc1 * x + Hc2) ** (-1 / 4.)


def u_exact(x):
    return Q0 / H_exact(x)

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import convert [as 别名]
def convert(value, u1, u2):
    "Convert value from units u1 to u2."
    unit_system = PISM.Context().unit_system
    return PISM.convert(unit_system, value, u1, u2)

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import convert [as 别名]
        vel_sia_observed.metadata(1).set_name('v_sia_observed')
        vel_sia_observed.metadata(1).set_string('long_name', "y-component of the 'observed' SIA velocities")

        vel_surface_observed = PISM.model.create2dVelocityVec(grid, "_surface_observed",
                                                              "observed surface velocities",
                                                              stencil_width=1)
        vel_surface_observed.copy_from(vel_sia_observed)
        vel_surface_observed.add(1., vel_ssa_observed)
        vecs.markForWriting(vel_surface_observed)
        final_velocity = vel_surface_observed

    sys = grid.ctx().unit_system()
    # Add the misfit weight.
    if misfit_weight_type == "fast":
        misfit_weight = fastIceMisfitWeight(modeldata, vel_ssa,
                                            PISM.convert(sys, fast_ice_speed, "m/year", "m/second"))
    else:
        misfit_weight = groundedIceMisfitWeight(modeldata)
    modeldata.vecs.add(misfit_weight, writing=True)

    if not noise is None:
        u_noise = PISM.vec.randVectorV(grid, noise / math.sqrt(2), final_velocity.get_stencil_width())
        final_velocity.add(PISM.convert(sys, 1.0, "m/year", "m/second"),
                           u_noise)

    pio = PISM.PIO(grid.com, "netcdf3")
    pio.open(output_file_name, 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(),