Python PISM.verbPrintf方法代码示例

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

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

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

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

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

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import verbPrintf [as 别名]
  def solve(self):
    if not self.ssa_init:
      pismVars = PISM.PISMVars()
      for var in [self.surface,self.thickness,self.bed,self.tauc,self.enthalpy,self.ice_mask]:
        pismVars.add(var)
        
      # The SSA instance will not keep a reference to pismVars; it only uses it to extract
      # its desired variables.  So it is safe to pass it pismVars and then let pismVars
      # go out of scope at the end of this method.
      self.ssa.init(pismVars)

      if not self.vel_bc is None:
        self.ssa.set_boundary_conditions(self.bc_mask,self.vel_bc)
      self.ssa_init = True
    
    PISM.verbPrintf(2,self.grid.com,"* Solving the SSA stress balance ...\n");
    fast = False;
    self.ssa.update(fast);
    self.vel_ssa = self.ssa.get_advective_2d_velocity()

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import verbPrintf [as 别名]
 def __call__(self, inverse_solver, count, data):
     """
     :param inverse_sovler: the solver (e.g. :class:`~InvSolver_Tikhonov`) we are listening to.
     :param count: the iteration number.
     :param data: dictionary of data related to the iteration.
     """
     method = inverse_solver.method
     if method != 'sd' and method != 'nlcg' and method != 'ign':
         if not self.didWarning:
             PISM.verbPrintf(1, PISM.Context().com, '\nWarning: unable to monitor adjoint for inverse method: %s\nOption -inv_monitor_adjoint ignored\n' % method)
         self.didWarning = True
         return
     fp = inverse_solver.forward_problem
     d = PISM.invert.sipletools.PISMLocalVector(data.d)
     r = PISM.invert.sipletools.PISMLocalVector(data.r)
     self.Td = fp.T(d, self.Td)
     self.TStarR = fp.TStar(r, out=self.TStarR)
     ip1 = fp.domainIP(d, self.TStarR)
     ip2 = fp.rangeIP(self.Td, r)
     logMessage("adjoint test: <Td,r>=%g <d,T^*r>=%g (percent error %g)", ip1, ip2, (abs(ip1 - ip2)) / max(abs(ip1), abs(ip2)))

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import verbPrintf [as 别名]
def run():
    context = PISM.Context()
    config = context.config
    com = context.com
    PISM.set_abort_on_sigint(True)

    WIDE_STENCIL = int(config.get_double("grid_max_stencil_width"))

    usage = \
        """  pismi.py [-i IN.nc [-o OUT.nc]]/[-a INOUT.nc] [-inv_data inv_data.nc] [-inv_forward model] 
                [-inv_design design_var] [-inv_method meth] 
    where:
    -i            IN.nc       is input file in NetCDF format: contains PISM-written model state
    -o            OUT.nc      is output file in NetCDF format to be overwritten
    -a            INOUT.nc    is input/output file in NetCDF format to be appended to
    -inv_data     inv_data.nc is data file containing extra inversion data (e.g. observed surface velocities)
    -inv_forward  model       forward model: only 'ssa' supported
    -inv_design   design_var  design variable name; one of 'tauc'/'hardav' for SSA inversions
    -inv_method   meth        algorithm for inversion [sd,nlcg,ign,tikhonov_lmvm]

    notes:
      * only one of -i/-a is allowed; both specify the input file
      * only one of -o/-a is allowed; both specify the output file
      * if -o is used, only the variables involved in inversion are written to the output file.
      * if -a is used, the varaibles involved in inversion are appended to the given file. No
        original variables in the file are changed.
   """

    append_mode = False
    PISM.setVerbosityLevel(1)

    input_filename = PISM.optionsString("-i", "input file")
    append_filename = PISM.optionsString("-a", "append file", default=None)
    output_filename = PISM.optionsString("-o", "output file", default=None)

    if (input_filename is None) and (append_filename is None):
        PISM.verbPrintf(1, com, "\nError: No input file specified. Use one of -i [file.nc] or -a [file.nc].\n")
        sys.exit(0)

    if (input_filename is not None) and (append_filename is not None):
        PISM.verbPrintf(1, com, "\nError: Only one of -i/-a is allowed.\n")
        sys.exit(0)

    if (output_filename is not None) and (append_filename is not None):
        PISM.verbPrintf(1, com, "\nError: Only one of -a/-o is allowed.\n")
        sys.edit(0)

    if append_filename is not None:
        input_filename = append_filename
        output_filename = append_filename
        append_mode = True

    inv_data_filename = PISM.optionsString("-inv_data", "inverse data file", default=input_filename)
    verbosity = PISM.optionsInt("-verbose", "verbosity level", default=2)

    do_plotting = PISM.optionsFlag("-inv_plot", "perform visualization during the computation", default=False)
    do_final_plot = PISM.optionsFlag("-inv_final_plot", "perform visualization at the end of the computation", default=False)
    Vmax = PISM.optionsReal("-inv_plot_vmax", "maximum velocity for plotting residuals", default=30)

    design_var = PISM.optionsList("-inv_ssa",
                                  "design variable for inversion",
                                  "tauc,hardav", "tauc")
    do_pause = PISM.optionsFlag("-inv_pause", "pause each iteration", default=False)

    do_restart = PISM.optionsFlag("-inv_restart", "Restart a stopped computation.", default=False)
    use_design_prior = PISM.optionsFlag("-inv_use_design_prior", "Use prior from inverse data file as initial guess.", default=True)

    prep_module = PISM.optionsString("-inv_prep_module", "Python module used to do final setup of inverse solver", default=None)

    is_regional = PISM.optionsFlag("-regional", "Compute SIA/SSA using regional model semantics", default=False)

    using_zeta_fixed_mask = PISM.optionsFlag("-inv_use_zeta_fixed_mask",
                                             "Enforce locations where the parameterized design variable should be fixed. (Automatically determined if not provided)", default=True)

    inv_method = config.get_string("inv_ssa_method")

    if output_filename is None:
        output_filename = "pismi_" + os.path.basename(input_filename)

    saving_inv_data = (inv_data_filename != output_filename)

    PISM.setVerbosityLevel(verbosity)
    forward_run = SSAForwardRun(input_filename, inv_data_filename, design_var)
    forward_run.setup()
    design_param = forward_run.designVariableParameterization()
    solver = PISM.invert.ssa.createInvSSASolver(forward_run)

    modeldata = forward_run.modeldata
    vecs = modeldata.vecs
    grid = modeldata.grid

    # Determine the prior guess for tauc/hardav. This can be one of
    # a) tauc/hardav from the input file (default)
    # b) tauc/hardav_prior from the inv_datafile if -inv_use_design_prior is set
    design_prior = createDesignVec(grid, design_var, '%s_prior' % design_var)
    long_name = design_prior.metadata().get_string("long_name")
    units = design_prior.metadata().get_string("units")
    design_prior.set_attrs("", "best prior estimate for %s (used for inversion)" % long_name, units, "")
    if PISM.util.fileHasVariable(inv_data_filename, "%s_prior" % design_var) and use_design_prior:
        PISM.logging.logMessage("  Reading '%s_prior' from inverse data file %s.\n" % (design_var, inv_data_filename))
#.........这里部分代码省略.........

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import verbPrintf [as 别名]
def test_lin(ssarun):
    grid = ssarun.grid

    PISM.verbPrintf(1, grid.com, "\nTest Linearization (Comparison with finite differences):\n")

    S = 250
    d_viewer = PETSc.Viewer().createDraw(title="d", size=S)
    Td_viewer = PETSc.Viewer().createDraw(title="Td", size=S)
    Td_fd_viewer = PETSc.Viewer().createDraw(title="Td_fd", size=S)
    d_Td_viewer = PETSc.Viewer().createDraw(title="d_Td", size=S)

    for (i, j) in grid.points():
        d = PISM.IceModelVec2S()
        d.create(grid, "", PISM.WITH_GHOSTS)
        d.set(0)
        with PISM.vec.Access(comm=d):
            d[i, j] = 1

        ssarun.ssa.linearize_at(zeta1)
        u1 = PISM.IceModelVec2V()
        u1.create(grid, "", PISM.WITH_GHOSTS)
        u1.copy_from(ssarun.ssa.solution())

        Td = PISM.IceModelVec2V()
        Td.create(grid, "", PISM.WITH_GHOSTS)
        ssarun.ssa.apply_linearization(d, Td)

        eps = 1e-8
        zeta2 = PISM.IceModelVec2S()
        zeta2.create(grid, "", PISM.WITH_GHOSTS)
        zeta2.copy_from(d)
        zeta2.scale(eps)
        zeta2.add(1, zeta1)
        ssarun.ssa.linearize_at(zeta2)
        u2 = PISM.IceModelVec2V()
        u2.create(grid, "", PISM.WITH_GHOSTS)
        u2.copy_from(ssarun.ssa.solution())

        Td_fd = PISM.IceModelVec2V()
        Td_fd.create(grid, "", PISM.WITH_GHOSTS)
        Td_fd.copy_from(u2)
        Td_fd.add(-1, u1)
        Td_fd.scale(1. / eps)

        d_Td = PISM.IceModelVec2V()
        d_Td.create(grid, "", PISM.WITH_GHOSTS)
        d_Td.copy_from(Td_fd)
        d_Td.add(-1, Td)

        n_Td_fd = Td_fd.norm(PETSc.NormType.NORM_2)
        n_Td_l2 = Td.norm(PETSc.NormType.NORM_2)
        n_Td_l1 = Td.norm(PETSc.NormType.NORM_1)
        n_Td_linf = Td.norm(PETSc.NormType.NORM_INFINITY)

        n_d_Td_l2 = d_Td.norm(PETSc.NormType.NORM_2)
        n_d_Td_l1 = d_Td.norm(PETSc.NormType.NORM_1)
        n_d_Td_linf = d_Td.norm(PETSc.NormType.NORM_INFINITY)

        PISM.verbPrintf(1, grid.com, "(i,j)=(%d,%d)\n" % (i, j))
        PISM.verbPrintf(1, grid.com, "apply_linearization(d): l2 norm %.10g; finite difference %.10g\n" % (n_Td_l2, n_Td_fd))

        r_d_l2 = 0
        if n_Td_l2 != 0:
            r_d_l2 = n_d_Td_l2 / n_Td_l2
        r_d_l1 = 0
        if n_Td_l1 != 0:
            r_d_l1 = n_d_Td_l1 / n_Td_l1

        r_d_linf = 0
        if n_Td_linf != 0:
            r_d_linf = n_d_Td_linf / n_Td_linf
        PISM.verbPrintf(1, grid.com, "relative difference: l2 norm %.10g l1 norm %.10g linf norm %.10g\n" % (r_d_l2, r_d_l1, r_d_linf))

        PISM.verbPrintf(1, grid.com, "\n")

        d_global = PISM.IceModelVec2S()
        d_global.create(grid, "", PISM.WITHOUT_GHOSTS)
        d_global.copy_from(d)
        d_global.get_vec().view(d_viewer)

        Td_global = PISM.IceModelVec2V()
        Td_global.create(grid, "", PISM.WITHOUT_GHOSTS)
        Td_global.copy_from(Td)
        # if n_Td_linf > 0:
        #   Td_global.scale(1./n_Td_linf)
        Td_global.get_vec().view(Td_viewer)

        Td_fd_global = PISM.IceModelVec2V()
        Td_fd_global.create(grid, "", PISM.WITHOUT_GHOSTS)
        Td_fd_global.copy_from(Td_fd)
        # if n_Td_fd_linf > 0:
        #   Td_fd_global.scale(1./n_Td_linf)
        Td_fd_global.get_vec().view(Td_fd_viewer)

        d_Td_global = PISM.IceModelVec2V()
        d_Td_global.create(grid, "", PISM.WITHOUT_GHOSTS)
        d_Td_global.copy_from(d_Td)
        # if n_Td_linf > 0:
        #   d_Td_global.scale(1./n_Td_linf)
        d_Td_global.get_vec().view(d_Td_viewer)
#.........这里部分代码省略.........

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import verbPrintf [as 别名]
def test_linearization_transpose(ssarun):
    grid = ssarun.grid

    S = 250
    r_viewer = PETSc.Viewer().createDraw(title="r", size=S)
    TStarR_viewer = PETSc.Viewer().createDraw(title="TStarR", size=S)
    TStarR_indirect_viewer = PETSc.Viewer().createDraw(title="TStarR (ind)", size=S)
    d_TStarR_viewer = PETSc.Viewer().createDraw(title="d_TStarR_fd", size=S)

    ssarun.ssa.linearize_at(zeta1)
    u = PISM.IceModelVec2V()
    u.create(grid, "", PISM.WITH_GHOSTS)
    u.copy_from(ssarun.ssa.solution())

    Td = PISM.IceModelVec2V()
    Td.create(grid, "", PISM.WITHOUT_GHOSTS)

    TStarR = PISM.IceModelVec2S()
    TStarR.create(grid, "", PISM.WITHOUT_GHOSTS)

    TStarR_indirect = PISM.IceModelVec2S()
    TStarR_indirect.create(grid, "", PISM.WITHOUT_GHOSTS)

    for (i, j) in grid.points():

        for k in range(2):

            r = PISM.IceModelVec2V()
            r.create(grid, "", PISM.WITH_GHOSTS)
            r.set(0)
            with PISM.vec.Access(comm=r):
                if k == 0:
                    r[i, j].u = 1
                else:
                    r[i, j].v = 1

            ssarun.ssa.apply_linearization_transpose(r, TStarR)

            r_global = PISM.IceModelVec2V()
            r_global.create(grid, "", PISM.WITHOUT_GHOSTS)
            r_global.copy_from(r)

            for (k, l) in grid.points():
                with PISM.vec.Access(nocomm=TStarR_indirect):
                    d = PISM.IceModelVec2S()
                    d.create(grid, "", PISM.WITH_GHOSTS)
                    d.set(0)
                    with PISM.vec.Access(comm=d):
                        d[k, l] = 1

                    ssarun.ssa.apply_linearization(d, Td)

                    TStarR_indirect[k, l] = Td.get_vec().dot(r_global.get_vec())

            d_TStarR = PISM.IceModelVec2S()
            d_TStarR.create(grid, "", PISM.WITHOUT_GHOSTS)

            d_TStarR.copy_from(TStarR)
            d_TStarR.add(-1, TStarR_indirect)

            PISM.verbPrintf(1, grid.com, "\nTest Linearization Transpose (%d,%d):\n" % (i, j))

            view(r_global, r_viewer)
            view(TStarR, TStarR_viewer)
            view(TStarR_indirect, TStarR_indirect_viewer)
            view(d_TStarR, d_TStarR_viewer)

            PISM.logging.pause()

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import verbPrintf [as 别名]
def test_j_design(ssarun):
    grid = ssarun.grid

    S = 250
    d_viewer = PETSc.Viewer().createDraw(title="d", size=S)
    drhs_viewer = PETSc.Viewer().createDraw(title="drhs", size=S)
    drhs_fd_viewer = PETSc.Viewer().createDraw(title="drhs_fd", size=S)
    d_drhs_viewer = PETSc.Viewer().createDraw(title="d_drhs", size=S)

    ssarun.ssa.linearize_at(zeta1)
    u1 = PISM.IceModelVec2V()
    u1.create(grid, "", PISM.WITH_GHOSTS)
    u1.copy_from(ssarun.ssa.solution())

    for (i, j) in grid.points():
        d = PISM.IceModelVec2S()
        d.create(grid, "", PISM.WITH_GHOSTS)
        d.set(0)
        with PISM.vec.Access(comm=d):
            d[i, j] = 1

        ssarun.ssa.linearize_at(zeta1)

        rhs1 = PISM.IceModelVec2V()
        rhs1.create(grid, "", PISM.WITHOUT_GHOSTS)
        ssarun.ssa.assemble_residual(u1, rhs1)

        eps = 1e-8
        zeta2 = PISM.IceModelVec2S()
        zeta2.create(grid, "zeta_prior", PISM.WITH_GHOSTS)
        zeta2.copy_from(d)
        zeta2.scale(eps)
        zeta2.add(1, zeta1)
        ssarun.ssa.set_design(zeta2)

        rhs2 = PISM.IceModelVec2V()
        rhs2.create(grid, "", PISM.WITHOUT_GHOSTS)
        ssarun.ssa.assemble_residual(u1, rhs2)

        drhs_fd = PISM.IceModelVec2V()
        drhs_fd.create(grid, "", PISM.WITHOUT_GHOSTS)
        drhs_fd.copy_from(rhs2)
        drhs_fd.add(-1, rhs1)
        drhs_fd.scale(1. / eps)

        drhs = PISM.IceModelVec2V()
        drhs.create(grid, "", PISM.WITHOUT_GHOSTS)
        ssarun.ssa.apply_jacobian_design(u1, d, drhs)

        d_drhs = PISM.IceModelVec2V()
        d_drhs.create(grid, "", PISM.WITHOUT_GHOSTS)

        d_drhs.copy_from(drhs)
        d_drhs.add(-1, drhs_fd)

        n_drhs_fd = drhs_fd.norm(PETSc.NormType.NORM_2)
        n_drhs_l2 = drhs.norm(PETSc.NormType.NORM_2)
        n_drhs_l1 = drhs.norm(PETSc.NormType.NORM_1)
        n_drhs_linf = drhs.norm(PETSc.NormType.NORM_INFINITY)

        n_d_drhs_l2 = d_drhs.norm(PETSc.NormType.NORM_2)
        n_d_drhs_l1 = d_drhs.norm(PETSc.NormType.NORM_1)
        n_d_drhs_linf = d_drhs.norm(PETSc.NormType.NORM_INFINITY)

        PISM.verbPrintf(1, grid.com, "\nTest Jacobian Design (Comparison with finite differences):\n")
        PISM.verbPrintf(1, grid.com, "jacobian_design(d): l2 norm %.10g; finite difference %.10g\n" % (n_drhs_l2, n_drhs_fd))
        if n_drhs_linf == 0:
            PISM.verbPrintf(1, grid.com, "difference: l2 norm %.10g l1 norm %.10g linf norm %.10g\n" % (n_d_drhs_l2, n_d_drhs_l1, n_d_drhs_linf))
        else:
            PISM.verbPrintf(1, grid.com, "relative difference: l2 norm %.10g l1 norm %.10g linf norm %.10g\n" % (n_d_drhs_l2 / n_drhs_l2, n_d_drhs_l1 / n_drhs_l1, n_d_drhs_linf / n_drhs_linf))

        view(d, d_viewer)
        view(drhs, drhs_viewer)
        view(drhs_fd, drhs_fd_viewer)
        view(d_drhs, d_drhs_viewer)

        PISM.logging.pause()

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import verbPrintf [as 别名]
    vecs = ssarun.modeldata.vecs
    grid = ssarun.grid

    # Determine the prior guess for tauc/hardav. This can be one of
    # a) tauc/hardav from the input file (default)
    # b) tauc/hardav_prior from the inv_datafile if -inv_use_design_prior is set
    design_prior = createDesignVec(grid, design_var, '%s_prior' % design_var)
    long_name = design_prior.string_attr("long_name")
    units = design_prior.string_attr("units")
    design_prior.set_attrs("", "best prior estimate for %s (used for inversion)" % long_name, units, "")
    if PISM.util.fileHasVariable(inv_data_filename, "%s_prior" % design_var) and use_design_prior:
        PISM.logging.logMessage("  Reading '%s_prior' from inverse data file %s.\n" % (design_var, inv_data_filename))
        design_prior.regrid(inv_data_filename, critical=True)
    else:
        if not PISM.util.fileHasVariable(input_filename, design_var):
            PISM.verbPrintf(1, com, "Initial guess for design variable is not available as '%s' in %s.\nYou can provide an initial guess in the inverse data file.\n" % (design_var, input_filename))
            exit(1)
        PISM.logging.logMessage("Reading '%s_prior' from '%s' in input file.\n" % (design_var, design_var))
        design = createDesignVec(grid, design_var)
        design.regrid(input_filename, True)
        design_prior.copy_from(design)

    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)

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import verbPrintf [as 别名]
  def report(self):
      grid = self.grid

      ssa_stdout = self.solver.ssa.stdout_report()
      PISM.verbPrintf(3,grid.com,ssa_stdout)

      maxvecerr = 0.0; avvecerr = 0.0; 
      avuerr = 0.0; avverr = 0.0;
      maxuerr = 0.0; maxverr = 0.0;

      if(self.config.get_flag("do_pseudo_plastic_till")):
        PISM.verbPrintf(1,grid.com, "WARNING: numerical errors not valid for pseudo-plastic till\n")
      PISM.verbPrintf(1,grid.com, "NUMERICAL ERRORS in velocity relative to exact solution:\n")

      vel_ssa = self.solver.ssa.get_advective_2d_velocity()
      
      vel_ssa.begin_access()

      exactvelmax = 0; gexactvelmax = 0;
      for (i,j) in self.grid.points():
        x=grid.x[i]; y=grid.y[j]
        (uexact,vexact) = self.exactSolution(i,j,x,y);
        exactnormsq=math.sqrt(uexact*uexact+vexact*vexact);
        exactvelmax = max(exactnormsq,exactvelmax);
        solution = vel_ssa[i,j]
        uerr = abs(solution.u-uexact)
        verr = abs(solution.v-vexact)
        avuerr += uerr;
        avverr += verr;
        maxuerr = max(maxuerr,uerr);
        maxverr = max(maxverr,verr)
        vecerr = math.sqrt(uerr * uerr + verr * verr);
        maxvecerr = max(maxvecerr,vecerr);
        avvecerr = avvecerr + vecerr;

      vel_ssa.end_access();
      
      gexactvelmax = PISM.globalMax(exactvelmax,grid.com);      
      gmaxuerr     = PISM.globalMax(maxuerr,grid.com);
      gmaxverr     = PISM.globalMax(maxverr,grid.com);
      gavuerr      = PISM.globalSum(avuerr,grid.com) / (grid.Mx*grid.My)
      gavverr      = PISM.globalSum(avverr,grid.com) / (grid.Mx*grid.My)
      gmaxvecerr   = PISM.globalMax(maxvecerr,grid.com)
      gavvecerr    = PISM.globalMax(avvecerr,grid.com) / (grid.Mx*grid.My)

      report_velocity_scale = PISM.secpera
      PISM.verbPrintf(1,grid.com, "velocity  :  maxvector   prcntavvec      maxu      maxv       avu       avv\n");
      #FIXME: variable arguments to verbPrintf are not working.  For now, do the string formatting on the python side.  Maybe
      #this is the best approach.
      PISM.verbPrintf(1,grid.com, "           %11.4f%13.5f%10.4f%10.4f%10.4f%10.4f\n",
                      gmaxvecerr*report_velocity_scale, (gavvecerr/gexactvelmax)*100.0,
                      gmaxuerr*report_velocity_scale, gmaxverr*report_velocity_scale, gavuerr*report_velocity_scale, 
                      gavverr*report_velocity_scale) 
      PISM.verbPrintf(1,grid.com, "NUM ERRORS DONE\n")

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import verbPrintf [as 别名]
def printing_test():
    "Test verbPrintf"
    ctx = PISM.Context()
    PISM.verbPrintf(1, ctx.com, "hello %s!\n", "world")

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import verbPrintf [as 别名]
def print_logger(message, verbosity):
    """Implements a logger that forwards messages to :cpp:func:`verbPrintf`."""
    com = PISM.Context().com
    msg = str(message)
    PISM.verbPrintf(verbosity, com, msg)