tIPFunctional.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
  HTML git clone git://src.adamsgaard.dk/pism
   DIR Log
   DIR Files
   DIR Refs
   DIR LICENSE
       ---
       tIPFunctional.cc (2818B)
       ---
            1 // Copyright (C) 2012, 2013, 2014, 2015, 2017  David Maxwell
            2 //
            3 // This file is part of PISM.
            4 //
            5 // PISM is free software; you can redistribute it and/or modify it under the
            6 // terms of the GNU General Public License as published by the Free Software
            7 // Foundation; either version 3 of the License, or (at your option) any later
            8 // version.
            9 //
           10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
           11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
           12 // FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
           13 // details.
           14 //
           15 // You should have received a copy of the GNU General Public License
           16 // along with PISM; if not, write to the Free Software
           17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
           18 
           19 #include "IPFunctional.hh"
           20 #include "pism/util/IceGrid.hh"
           21 #include "pism/util/error_handling.hh"
           22 
           23 namespace pism {
           24 namespace inverse {
           25 
           26 void gradientFD(IPFunctional<IceModelVec2S> &f, IceModelVec2S &x, IceModelVec2S &gradient) {
           27   const IceGrid &grid = *x.grid();
           28   double h = PETSC_SQRT_MACHINE_EPSILON;
           29 
           30   double F0,Fh;
           31   
           32   f.valueAt(x,&F0);
           33   
           34   IceModelVec::AccessList list(gradient);
           35 
           36   ParallelSection loop(grid.com);
           37   try {
           38     for (Points p(grid); p; p.next()) {
           39       const int i = p.i(), j = p.j();
           40 
           41       {
           42         IceModelVec::AccessList access_x(x);
           43         x(i,j) += h;
           44       }
           45       x.update_ghosts();
           46 
           47       f.valueAt(x,&Fh);
           48 
           49       {
           50         IceModelVec::AccessList access_x(x);
           51         x(i,j) -= h;
           52       }
           53       x.update_ghosts();
           54 
           55       gradient(i,j) = (Fh-F0)/h;
           56     }
           57   } catch (...) {
           58     loop.failed();
           59   }
           60   loop.check();
           61 }
           62 
           63 void gradientFD(IPFunctional<IceModelVec2V> &f, IceModelVec2V &x, IceModelVec2V &gradient) {
           64   const IceGrid &grid = *x.grid();
           65   double h = PETSC_SQRT_MACHINE_EPSILON;
           66 
           67   double F0,Fh;
           68   
           69   f.valueAt(x,&F0);
           70   
           71   IceModelVec::AccessList access_gradient(gradient);
           72 
           73   ParallelSection loop(grid.com);
           74   try {
           75     for (Points p(grid); p; p.next()) {
           76       const int i = p.i(), j = p.j();
           77 
           78       {
           79         IceModelVec::AccessList access_x(x);
           80         x(i,j).u += h;
           81       }
           82       x.update_ghosts();
           83 
           84       f.valueAt(x,&Fh);
           85 
           86       {
           87         IceModelVec::AccessList access_x(x);
           88         x(i,j).u -= h;
           89       }
           90       x.update_ghosts();
           91 
           92       gradient(i,j).u = (Fh-F0)/h;
           93 
           94       {
           95         IceModelVec::AccessList access_x(x);
           96         x(i,j).v += h;
           97       }
           98       x.update_ghosts();
           99 
          100       f.valueAt(x,&Fh);
          101 
          102       {
          103         IceModelVec::AccessList access_x(x);
          104         x(i,j).v -= h;
          105       }
          106       x.update_ghosts();
          107 
          108       gradient(i,j).v = (Fh-F0)/h;
          109     }
          110   } catch (...) {
          111     loop.failed();
          112   }
          113   loop.check();
          114 }
          115 
          116 // PetscErrorCode gradientFD(IPFunctional<IceModelVec2V> &f, IceModelVec2V &x, IceModelVec2V &gradient);
          117 
          118 } // end of namespace inverse
          119 } // end of namespace pism