Python PISM.globalSum方法代码示例

# 需要导入模块: import PISM [as 别名]
# 或者: from PISM import globalSum [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 globalSum [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 globalSum [as 别名]
  def report(self):
    uerr = 0.0; verr=0.0; relvecerr=0.0; accN=0.0 
    accArea=0.0; maxcComputed=0.0; vecErrAcc = 0.0;
    
    grid = self.grid
    area = grid.dx*grid.dy
  
    mask = self.solver.ice_mask;
    H = self.solver.thickness;
    azi=self.obsAzimuth;
    mag = self.obsMagnitude;
    acc = self.obsAccurate;
    vel_ssa = self.solver.solution();

    m = PISM.MaskQuery(mask)
    
    with PISM.util.Access([mask,H,azi,mag,acc,vel_ssa]):
      for (i,j) in grid.points():
        if m.ocean(i,j) and H[i,j] > 1.0:
          ccomputed = vel_ssa[i,j].magnitude()
          maxcComputed = max(maxcComputed,ccomputed)
          if( abs(acc[i,j]-1.0) < 0.1):
            accN += 1.0
            accArea += area
            uobs = mag[i,j] * math.sin((math.pi/180.0) * azi[i,j]);
            vobs = mag[i,j] * math.cos((math.pi/180.0) * azi[i,j]);
            Dv = abs(vobs-vel_ssa[i,j].v)
            Du = abs(uobs-vel_ssa[i,j].u)
            verr += Dv; uerr += Du
            relvecerr += (Dv*Dv+Du*Du) / (vobs*vobs+uobs*uobs)
            vecErrAcc += (Dv*Dv+Du*Du) * area

    gmaxcComputed = PISM.globalMax(maxcComputed,grid.com)
    gaccN         = PISM.globalSum(accN, grid.com)
    gaccArea      = PISM.globalSum(accArea, grid.com)
    gverr         = PISM.globalSum(verr, grid.com)
    guerr         = PISM.globalSum(uerr, grid.com)
    grelvecerr    = PISM.globalSum(relvecerr,grid.com)
    gvecErrAcc    = PISM.globalSum(vecErrAcc, grid.com)

    secpera = PISM.secpera
    r=PISM.VerbPrintf(self.grid.com,verbosity=2)
    r.println("maximum computed speed in ice shelf is %10.3f (m/a)", gmaxcComputed * secpera);
    r.println("ERRORS relative to observations of Ross Ice Shelf:");
    r.println("  [number of grid points in 'accurate observed area' = %d]", int(gaccN))
    r.println("  [area of 'accurate observed area' = %9.4f (km^2)]", gaccArea / 1e6)
    r.println("  following are average errors computed over 'accurate observed area':");
    r.println("  average error in x-comp of vel       = %9.3f (m/a)", (gverr * secpera) / gaccN);
    r.println("  average error in y-comp of vel       = %9.3f (m/a)",(guerr * secpera) / gaccN);
    r.println("  average relative error in vector vel = %9.5f", grelvecerr / gaccN)
    gvecErrAcc = secpera * math.sqrt(gvecErrAcc) / math.sqrt(gaccArea);
    r.println("  rms average error in vector vel      = %9.3f (m/a)\n", gvecErrAcc);
    
    if not self.riggs_file is None:
      self.report_riggs()