stable0.3/src/base/iceModelVec.cc Source File

00001 // Copyright (C) 2008--2010 Ed Bueler and Constantine Khroulev
00002 //
00003 // This file is part of PISM.
00004 //
00005 // PISM is free software; you can redistribute it and/or modify it under the
00006 // terms of the GNU General Public License as published by the Free Software
00007 // Foundation; either version 2 of the License, or (at your option) any later
00008 // version.
00009 //
00010 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
00011 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00012 // FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
00013 // details.
00014 //
00015 // You should have received a copy of the GNU General Public License
00016 // along with PISM; if not, write to the Free Software
00017 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00018 
00019 #include "pism_const.hh"
00020 #include "iceModelVec.hh"
00021 
00022 
00023 IceModelVec::IceModelVec() {
00024 
00025   shallow_copy = false;
00026   v = PETSC_NULL;
00027   da = PETSC_NULL;
00028   grid = PETSC_NULL;
00029   array = PETSC_NULL;
00030   localp = true;
00031   use_interpolation_mask = false;
00032 
00033   dims = GRID_2D;               // default
00034   dof = 1;                      // default
00035   s_width = 1;
00036 
00037   access_counter = 0;
00038 
00039   name = "unintialized variable";
00040 
00041   vars.resize(dof);
00042 
00043   map_viewers = new map<string,PetscViewer>;
00044   reset_attrs();
00045 }
00046 
00047 
00049 
00055 IceModelVec::IceModelVec(const IceModelVec &other) {
00056   // This IceModelVec is a shallow copy!
00057   shallow_copy = true;
00058 
00059   v = other.v;
00060   da = other.da;
00061   dims = other.dims;
00062   dof = other.dof;
00063   grid = other.grid;
00064   array = other.array;
00065   access_counter = other.access_counter;
00066   localp = other.localp;
00067 
00068   map_viewers = other.map_viewers;
00069 
00070   use_interpolation_mask = other.use_interpolation_mask;
00071   interpolation_mask = other.interpolation_mask;
00072 
00073   name = other.name;
00074 
00075   vars = other.vars;
00076 }
00077 
00078 
00079 IceModelVec::~IceModelVec() {
00080   // Only destroy the IceModelVec if it is not a shallow copy:
00081   if (!shallow_copy) destroy();
00082 }
00083 
00085 bool IceModelVec::was_created() {
00086   return (v != PETSC_NULL);
00087 }
00088 
00090 GridType IceModelVec::grid_type() {
00091   return dims;
00092 }
00093 
00094 PetscErrorCode  IceModelVec::destroy() {
00095   PetscErrorCode ierr;
00096 
00097   if (v != PETSC_NULL) {
00098     ierr = VecDestroy(v); CHKERRQ(ierr);
00099     v = PETSC_NULL;
00100   }
00101   if (da != PETSC_NULL) {
00102     ierr = DADestroy(da); CHKERRQ(ierr);
00103     da = PETSC_NULL;
00104   }
00105 
00106   // map-plane viewers:
00107   if (map_viewers != NULL) {
00108     map<string,PetscViewer>::iterator i;
00109     for (i = (*map_viewers).begin(); i != (*map_viewers).end(); ++i) {
00110       if ((*i).second != PETSC_NULL) {
00111         ierr = PetscViewerDestroy((*i).second); CHKERRQ(ierr);
00112       }
00113     }
00114     delete map_viewers;
00115     map_viewers = NULL;
00116   }
00117 
00118   return 0;
00119 }
00120 
00122 
00128 PetscErrorCode IceModelVec::range(PetscReal &min, PetscReal &max) {
00129   PetscReal my_min, my_max, gmin, gmax;
00130   PetscErrorCode ierr;
00131   ierr = checkAllocated(); CHKERRQ(ierr);
00132 
00133   ierr = VecMin(v, PETSC_NULL, &my_min); CHKERRQ(ierr);
00134   ierr = VecMax(v, PETSC_NULL, &my_max); CHKERRQ(ierr);
00135 
00136   if (localp) {
00137     // needs a reduce operation; use PetscGlobalMax;
00138     ierr = PetscGlobalMin(&my_min, &gmin, grid->com); CHKERRQ(ierr);
00139     ierr = PetscGlobalMax(&my_max, &gmax, grid->com); CHKERRQ(ierr);
00140     min = gmin;
00141     max = gmax;
00142   } else {
00143     min = my_min;
00144     max = my_max;
00145   }
00146   return 0;
00147 }
00148 
00149 
00151 
00155 PetscErrorCode IceModelVec::norm(NormType n, PetscReal &out) {
00156   PetscErrorCode ierr;
00157   PetscReal my_norm, gnorm;
00158   ierr = checkAllocated(); CHKERRQ(ierr);
00159 
00160   ierr = VecNorm(v, n, &my_norm); CHKERRQ(ierr);
00161 
00162   if (localp) {
00163     // needs a reduce operation; use PetscGlobalMax if NORM_INFINITY,
00164     //   otherwise PetscGlobalSum; carefully in NORM_2 case
00165     if (n == NORM_1_AND_2) {
00166       SETERRQ1(1, 
00167          "IceModelVec::norm(...): NORM_1_AND_2 not implemented (called as %s.norm(...))\n",
00168          name.c_str());
00169     } else if (n == NORM_1) {
00170       ierr = PetscGlobalSum(&my_norm, &gnorm, grid->com); CHKERRQ(ierr);
00171     } else if (n == NORM_2) {
00172       my_norm = PetscSqr(my_norm);  // undo sqrt in VecNorm before sum
00173       ierr = PetscGlobalSum(&my_norm, &gnorm, grid->com); CHKERRQ(ierr);
00174       gnorm = sqrt(gnorm);
00175     } else if (n == NORM_INFINITY) {
00176       ierr = PetscGlobalMax(&my_norm, &gnorm, grid->com); CHKERRQ(ierr);
00177     } else {
00178       SETERRQ1(2, "IceModelVec::norm(...): unknown NormType (called as %s.norm(...))\n",
00179          name.c_str());
00180     }
00181     out = gnorm;
00182   } else {
00183     out = my_norm;
00184   }
00185   return 0;
00186 }
00187 
00188 
00190 
00193 PetscErrorCode IceModelVec::squareroot() {
00194   PetscErrorCode ierr;
00195   ierr = checkAllocated(); CHKERRQ(ierr);
00196 
00197   ierr = VecSqrt(v); CHKERRQ(ierr);
00198   return 0;
00199 }
00200 
00201 
00203 PetscErrorCode IceModelVec::add(PetscScalar alpha, IceModelVec &x) {
00204   PetscErrorCode ierr;
00205   ierr = checkAllocated(); CHKERRQ(ierr);
00206   ierr = x.checkAllocated(); CHKERRQ(ierr);
00207   ierr = checkCompatibility("add", x); CHKERRQ(ierr);
00208 
00209   ierr = VecAXPY(v, alpha, x.v); CHKERRQ(ierr);
00210   return 0;
00211 }
00212 
00214 PetscErrorCode IceModelVec::add(PetscScalar alpha, IceModelVec &x, IceModelVec &result) {
00215   PetscErrorCode ierr;
00216   ierr = checkAllocated(); CHKERRQ(ierr);
00217   ierr = x.checkAllocated(); CHKERRQ(ierr);
00218   ierr = result.checkAllocated(); CHKERRQ(ierr);
00219   ierr = checkCompatibility("add", x); CHKERRQ(ierr);
00220   ierr = checkCompatibility("add", result); CHKERRQ(ierr);
00221   
00222   ierr = VecWAXPY(result.v, alpha, x.v, v); CHKERRQ(ierr);
00223   return 0;
00224 }
00225 
00227 PetscErrorCode IceModelVec::shift(PetscScalar alpha) {
00228   PetscErrorCode ierr;
00229   ierr = checkAllocated(); CHKERRQ(ierr);
00230 
00231   ierr = VecShift(v, alpha); CHKERRQ(ierr);
00232   return 0;
00233 }
00234 
00236 PetscErrorCode IceModelVec::scale(PetscScalar alpha) {
00237   PetscErrorCode ierr;
00238   ierr = checkAllocated(); CHKERRQ(ierr);
00239 
00240   ierr = VecScale(v, alpha); CHKERRQ(ierr);
00241   return 0;
00242 }
00243 
00245 PetscErrorCode  IceModelVec::multiply_by(IceModelVec &x, IceModelVec &result) {
00246   PetscErrorCode ierr;
00247   ierr = checkAllocated(); CHKERRQ(ierr);
00248   ierr = x.checkAllocated(); CHKERRQ(ierr);
00249   ierr = checkCompatibility("multiply_by", x); CHKERRQ(ierr);
00250   ierr = checkCompatibility("multiply_by", result); CHKERRQ(ierr);
00251 
00252   ierr = VecPointwiseMult(result.v, v, x.v); CHKERRQ(ierr);
00253   return 0;
00254 }
00255 
00257 PetscErrorCode  IceModelVec::multiply_by(IceModelVec &x) {
00258   PetscErrorCode ierr;
00259   ierr = checkAllocated(); CHKERRQ(ierr);
00260   ierr = x.checkAllocated(); CHKERRQ(ierr);
00261   ierr = checkCompatibility("multiply_by", x); CHKERRQ(ierr);
00262 
00263   ierr = VecPointwiseMult(v, v, x.v); CHKERRQ(ierr);
00264   return 0;
00265 }
00266 
00268 
00271 PetscErrorCode  IceModelVec::copy_to(Vec destination) {
00272   PetscErrorCode ierr;
00273   ierr = checkAllocated(); CHKERRQ(ierr);
00274   if (localp) {
00275     ierr = DALocalToGlobal(da, v, INSERT_VALUES, destination); CHKERRQ(ierr);
00276   } else {
00277     ierr = VecCopy(v, destination); CHKERRQ(ierr);
00278   }
00279   return 0;
00280 }
00281 
00282 PetscErrorCode IceModelVec::copy_from(Vec source) {
00283   PetscErrorCode ierr;
00284   ierr = checkAllocated(); CHKERRQ(ierr);
00285   if (localp) {
00286     ierr =   DAGlobalToLocalBegin(da, source, INSERT_VALUES, v);  CHKERRQ(ierr);
00287     ierr =     DAGlobalToLocalEnd(da, source, INSERT_VALUES, v);  CHKERRQ(ierr);
00288   } else {
00289     ierr = VecCopy(source, v); CHKERRQ(ierr);
00290   }
00291   return 0;
00292 }
00293 
00294 
00296 PetscErrorCode  IceModelVec::copy_to(IceModelVec &destination) {
00297   PetscErrorCode ierr;
00298   ierr = checkAllocated(); CHKERRQ(ierr);
00299   ierr = destination.checkAllocated(); CHKERRQ(ierr);
00300   ierr = checkCompatibility("copy_to", destination); CHKERRQ(ierr);
00301 
00302   ierr = VecCopy(v, destination.v); CHKERRQ(ierr);
00303   return 0;
00304 }
00305 
00307 PetscErrorCode  IceModelVec::copy_from(IceModelVec &source) {
00308   PetscErrorCode ierr;
00309   ierr = checkAllocated(); CHKERRQ(ierr);
00310   ierr = source.checkAllocated(); CHKERRQ(ierr);
00311   ierr = checkCompatibility("copy_from", source); CHKERRQ(ierr);
00312 
00313   ierr = VecCopy(source.v, v); CHKERRQ(ierr);
00314   return 0;
00315 }
00316 
00318 PetscErrorCode  IceModelVec::set_name(const char new_name[], int N) {
00319   reset_attrs();
00320   
00321   if (N == 0)
00322     name = new_name;
00323 
00324   vars[N].short_name = new_name;
00325 
00326   return 0;
00327 }
00328 
00330 
00335 PetscErrorCode  IceModelVec::set_glaciological_units(string my_units) {
00336 
00337   PetscErrorCode ierr;
00338     
00339   for (int j = 0; j < dof; ++j) {
00340    ierr = vars[j].set_glaciological_units(my_units); CHKERRQ(ierr);
00341   }
00342   
00343   return 0;
00344 }
00345 
00347 PetscErrorCode IceModelVec::reset_attrs() {
00348 
00349   time_independent = false;
00350   write_in_glaciological_units = false;
00351   output_data_type = NC_DOUBLE;
00352 
00353   for (int j = 0; j < dof; ++j) {
00354     vars[j].reset();
00355     // _FillValue should be a number and should depend on the quantity in
00356     //     question. It should also be outside the valid range.
00357     //     vars[j].set("_FillValue", GSL_NAN);
00358   }
00359 
00360   return 0;
00361 }
00362 
00364 
00372 PetscErrorCode IceModelVec::set_attrs(string my_pism_intent,
00373                                       string my_long_name,
00374                                       string my_units,
00375                                       string my_standard_name,
00376                                       int N) {
00377 
00378   if (!my_long_name.empty()) {
00379     vars[N].set_string("long_name", my_long_name);
00380   }
00381 
00382   if (!my_units.empty()) {
00383     PetscErrorCode ierr = vars[N].set_units(my_units); CHKERRQ(ierr);
00384   }
00385 
00386   if (!my_pism_intent.empty()) {
00387     vars[N].set_string("pism_intent", my_pism_intent);
00388   }
00389 
00390   if (!my_standard_name.empty()) {
00391     vars[N].set_string("standard_name", my_standard_name);
00392   }
00393 
00394   return 0;
00395 }
00396 
00398 
00400 PetscErrorCode IceModelVec::regrid(const char filename[], LocalInterpCtx &lic, 
00401                                    bool critical) {
00402   PetscErrorCode ierr;
00403   MaskInterp *m = NULL;
00404   Vec g;
00405 
00406   if (dof != 1)
00407     SETERRQ(1, "This method only supports IceModelVecs with dof == 1.");
00408 
00409   if (use_interpolation_mask) m = &interpolation_mask;
00410 
00411   if (localp) {
00412     ierr = DACreateGlobalVector(da, &g); CHKERRQ(ierr);
00413 
00414     ierr = vars[0].regrid(filename, lic, critical, false, 0.0, m, g); CHKERRQ(ierr);
00415 
00416     ierr = DAGlobalToLocalBegin(da, g, INSERT_VALUES, v); CHKERRQ(ierr);
00417     ierr = DAGlobalToLocalEnd(da, g, INSERT_VALUES, v); CHKERRQ(ierr);
00418 
00419     ierr = VecDestroy(g); CHKERRQ(ierr);
00420   } else {
00421     ierr = vars[0].regrid(filename, lic, critical, false, 0.0, m, v); CHKERRQ(ierr);
00422   }
00423 
00424   return 0;
00425 }
00426 
00428 
00430 PetscErrorCode IceModelVec::regrid(const char filename[], 
00431                                    LocalInterpCtx &lic, PetscScalar default_value) {
00432   PetscErrorCode ierr;
00433   MaskInterp *m = NULL;
00434   Vec g;
00435 
00436   if (dof != 1)
00437     SETERRQ(1, "This method only supports IceModelVecs with dof == 1.");
00438 
00439   if (use_interpolation_mask) m = &interpolation_mask;
00440 
00441   if (localp) {
00442     ierr = DACreateGlobalVector(da, &g); CHKERRQ(ierr);
00443 
00444     ierr = vars[0].regrid(filename, lic, false, true, default_value, m, g); CHKERRQ(ierr);
00445 
00446     ierr = DAGlobalToLocalBegin(da, g, INSERT_VALUES, v); CHKERRQ(ierr);
00447     ierr = DAGlobalToLocalEnd(da, g, INSERT_VALUES, v); CHKERRQ(ierr);
00448 
00449     ierr = VecDestroy(g); CHKERRQ(ierr);
00450   } else {
00451     ierr = vars[0].regrid(filename, lic, false, true, default_value, m, v); CHKERRQ(ierr);
00452   }
00453 
00454   return 0;
00455 }
00456 
00458 PetscErrorCode IceModelVec::read(const char filename[], const unsigned int time) {           
00459   PetscErrorCode ierr;
00460   Vec g;
00461 
00462   if (getVerbosityLevel() > 3) {
00463     ierr = PetscPrintf(grid->com, "  Reading %s...\n", name.c_str()); CHKERRQ(ierr);
00464   }
00465 
00466   if (dof != 1)
00467     SETERRQ(1, "This method only supports IceModelVecs with dof == 1.");
00468 
00469   if (localp) {
00470     ierr = DACreateGlobalVector(da, &g); CHKERRQ(ierr);
00471 
00472     ierr = vars[0].read(filename, time, g); CHKERRQ(ierr);
00473 
00474     ierr = DAGlobalToLocalBegin(da, g, INSERT_VALUES, v); CHKERRQ(ierr);
00475     ierr = DAGlobalToLocalEnd(da, g, INSERT_VALUES, v); CHKERRQ(ierr);
00476 
00477     ierr = VecDestroy(g); CHKERRQ(ierr);
00478   } else {
00479     ierr = vars[0].read(filename, time, v); CHKERRQ(ierr);
00480   }
00481 
00482   return 0;
00483 }
00484 
00486 PetscErrorCode IceModelVec::write(const char filename[]) {
00487   PetscErrorCode ierr;
00488 
00489   if (getVerbosityLevel() > 3) {
00490     ierr = PetscPrintf(grid->com, "  Writing %s...\n", name.c_str()); CHKERRQ(ierr);
00491   }
00492   
00493   ierr = write(filename, output_data_type); CHKERRQ(ierr);
00494 
00495   return 0;
00496 }
00497 
00499 PetscErrorCode IceModelVec::write(const char filename[], nc_type nctype) {
00500   PetscErrorCode ierr;
00501   Vec g;
00502 
00503   if (dof != 1)
00504     SETERRQ(1, "This method only supports IceModelVecs with dof == 1.");
00505 
00506   vars[0].time_independent = time_independent;
00507 
00508   if (localp) {
00509     ierr = DACreateGlobalVector(da, &g); CHKERRQ(ierr);
00510     ierr = DALocalToGlobal(da, v, INSERT_VALUES, g); CHKERRQ(ierr);
00511 
00512     ierr = vars[0].write(filename, nctype, write_in_glaciological_units, g); CHKERRQ(ierr);
00513 
00514     ierr = VecDestroy(g); CHKERRQ(ierr);
00515   } else {
00516     ierr = vars[0].write(filename, nctype, write_in_glaciological_units, v); CHKERRQ(ierr);
00517   }
00518 
00519   return 0;
00520 }
00521 
00523 PetscErrorCode  IceModelVec::checkAllocated() {
00524   if (v == PETSC_NULL) {
00525     SETERRQ1(1,"IceModelVec ERROR: IceModelVec with name='%s' WAS NOT allocated\n",
00526              name.c_str());
00527   }
00528   return 0;
00529 }
00530 
00532 PetscErrorCode  IceModelVec::checkHaveArray() {
00533   PetscErrorCode ierr;
00534   ierr = checkAllocated(); CHKERRQ(ierr);
00535   if (array == PETSC_NULL) {
00536     SETERRQ1(1,"array for IceModelVec with name='%s' not available\n"
00537                "  (REMEMBER TO RUN begin_access() before access and end_access() after access)\n",
00538                name.c_str());
00539   }
00540   return 0;
00541 }
00542 
00544 PetscErrorCode IceModelVec::checkCompatibility(const char* func, IceModelVec &other) {
00545   PetscErrorCode ierr;
00546   PetscInt X_size, Y_size;
00547 
00548   if (dof != other.dof) {
00549     SETERRQ1(1, "IceModelVec::%s(...): operands have different numbers of degrees of freedom",
00550              func);
00551   }
00552 
00553   if (dims != other.dims) {
00554     SETERRQ1(1, "IceModelVec::%s(...): operands have different numbers of dimensions",
00555              func);
00556   }
00557 
00558   ierr = VecGetSize(v, &X_size); CHKERRQ(ierr);
00559   ierr = VecGetSize(other.v, &Y_size); CHKERRQ(ierr);
00560   if (X_size != Y_size)
00561     SETERRQ3(1, "IceModelVec::%s(...): incompatible Vec sizes (called as %s.%s(...))\n",
00562              func, name.c_str(), func);
00563 
00564 
00565   return 0;
00566 }
00567 
00569 PetscErrorCode  IceModelVec::begin_access() {
00570   PetscErrorCode ierr;
00571 #ifdef PISM_DEBUG
00572   ierr = checkAllocated(); CHKERRQ(ierr);
00573 
00574   if (access_counter < 0)
00575     SETERRQ(1, "IceModelVec::begin_access(): access_counter < 0");
00576 #endif
00577 
00578   if (access_counter == 0) {
00579     ierr = DAVecGetArray(da, v, &array); CHKERRQ(ierr);
00580   }
00581 
00582   access_counter++;
00583 
00584   return 0;
00585 }
00586 
00588 PetscErrorCode  IceModelVec::end_access() {
00589   PetscErrorCode ierr;
00590   access_counter--;
00591 
00592 #ifdef PISM_DEBUG
00593   ierr = checkAllocated(); CHKERRQ(ierr);
00594 
00595   if (access_counter < 0)
00596     SETERRQ(1, "IceModelVec::end_access(): access_counter < 0");
00597 #endif
00598 
00599   if (access_counter == 0) {
00600     ierr = DAVecRestoreArray(da, v, &array); CHKERRQ(ierr);
00601     array = NULL;
00602   }
00603 
00604   return 0;
00605 }
00606 
00608 PetscErrorCode  IceModelVec::beginGhostComm() {
00609   PetscErrorCode ierr;
00610   if (!localp) {
00611     SETERRQ1(1,"makes no sense to communicate ghosts for GLOBAL IceModelVec! (has name='%s')\n",
00612                name.c_str());
00613   }
00614   ierr = checkAllocated(); CHKERRQ(ierr);
00615   ierr = DALocalToLocalBegin(da, v, INSERT_VALUES, v);  CHKERRQ(ierr);
00616   return 0;
00617 }
00618 
00620 PetscErrorCode  IceModelVec::endGhostComm() {
00621   PetscErrorCode ierr;
00622   if (!localp) {
00623     SETERRQ1(1,"makes no sense to communicate ghosts for GLOBAL IceModelVec! (has name='%s')\n",
00624                name.c_str());
00625   }
00626   ierr = checkAllocated(); CHKERRQ(ierr);
00627   ierr = DALocalToLocalEnd(da, v, INSERT_VALUES, v); CHKERRQ(ierr);
00628   return 0;
00629 }
00630 
00632 PetscErrorCode  IceModelVec::beginGhostComm(IceModelVec &destination) {
00633   PetscErrorCode ierr;
00634   if (!localp) {
00635     SETERRQ1(1,"makes no sense to communicate ghosts for GLOBAL IceModelVec! (has name='%s')",
00636                name.c_str());
00637   }
00638 
00639   if (!destination.localp) {
00640     SETERRQ(1, "IceModelVec::beginGhostComm(): destination has to be local.");
00641   }
00642 
00643   if (dims != destination.dims) {
00644     SETERRQ(1, "IceModelVec::beginGhostComm(): operands have different numbers of dimensions.");
00645   }
00646 
00647   ierr = checkAllocated(); CHKERRQ(ierr);
00648   ierr = DALocalToLocalBegin(da, v, INSERT_VALUES, destination.v);  CHKERRQ(ierr);
00649   return 0;
00650 }
00651 
00653 PetscErrorCode  IceModelVec::endGhostComm(IceModelVec &destination) {
00654   PetscErrorCode ierr;
00655   if (!localp) {
00656     SETERRQ1(1,"makes no sense to communicate ghosts for GLOBAL IceModelVec! (has name='%s')\n",
00657                name.c_str());
00658   }
00659 
00660   if (!destination.localp) {
00661     SETERRQ(1, "IceModelVec::beginGhostComm(): destination has to be local.");
00662   }
00663 
00664   if (dims != destination.dims) {
00665     SETERRQ(1, "IceModelVec::beginGhostComm(): operands have different numbers of dimensions.");
00666   }
00667 
00668   ierr = checkAllocated(); CHKERRQ(ierr);
00669   ierr = DALocalToLocalEnd(da, v, INSERT_VALUES, destination.v); CHKERRQ(ierr);
00670   return 0;
00671 }
00672 
00674 PetscErrorCode  IceModelVec::set(const PetscScalar c) {
00675   PetscErrorCode ierr;
00676   ierr = checkAllocated(); CHKERRQ(ierr);
00677   ierr = VecSet(v,c); CHKERRQ(ierr);
00678   return 0;
00679 }
00680 
00682 
00685 bool IceModelVec::is_valid(PetscScalar a, int N) {
00686   return vars[N].is_valid(a);
00687 }
00688 
00690 
00692 PetscErrorCode IceModelVec::set_attr(string attr, string value, int N) {
00693   PetscErrorCode ierr;
00694 
00695   if (attr == "units") {
00696     ierr = vars[N].set_units(value); CHKERRQ(ierr);
00697     return 0;
00698   }
00699 
00700   if (attr == "glaciological_units") {
00701     ierr = vars[N].set_glaciological_units(value); CHKERRQ(ierr);
00702     return 0;
00703   }
00704 
00705   vars[N].set_string(attr, value);
00706   return 0;
00707 }
00708 
00710 PetscErrorCode IceModelVec::set_attr(string name, double value, int N) {
00711   vars[N].set(name, value);
00712   return 0;
00713 }
00714 
00716 PetscErrorCode IceModelVec::set_attr(string name, vector<double> values,
00717                                      int N) {
00718   vars[N].doubles[name] = values;
00719   return 0;
00720 }
00721 
00722 bool IceModelVec::has_attr(string name, int N) {
00723   return vars[N].has(name);
00724 }
00725 
00727 double IceModelVec::double_attr(string name, int N) {
00728   return vars[N].get(name);
00729 }
00730 
00732 string IceModelVec::string_attr(string n, int N) {
00733 
00734   if (n == "name")
00735     return name;
00736 
00737   return vars[N].get_string(n);
00738 }
00739 
00741 vector<double> IceModelVec::array_attr(string name, int N) {
00742   return vars[N].doubles[name];
00743 }
00744 
00746 PetscErrorCode IceModelVec::create_viewer(PetscInt viewer_size, string title, PetscViewer &viewer) {
00747   PetscErrorCode ierr;
00748   int x, y;
00749 
00750   ierr = compute_viewer_size(viewer_size, x, y); CHKERRQ(ierr);
00751 
00752   // note we reverse x <-> y; see IceGrid::createDA() for original reversal
00753   ierr = PetscViewerDrawOpen(grid->com, PETSC_NULL, title.c_str(),
00754                              PETSC_DECIDE, PETSC_DECIDE, y, x, &viewer);  CHKERRQ(ierr);
00755 
00756   // following should be redundant, but may put up a title even under 2.3.3-p1:3 where
00757   // there is a no titles bug
00758   PetscDraw draw;
00759   ierr = PetscViewerDrawGetDraw(viewer, 0, &draw); CHKERRQ(ierr);
00760   ierr = PetscDrawSetTitle(draw, title.c_str()); CHKERRQ(ierr);
00761 
00762   return 0;
00763 }
00764 
00766 PetscErrorCode IceModelVec::compute_viewer_size(int target_size, int &x, int &y) {
00767   PetscScalar Lx = grid->Lx, Ly = grid->Ly;
00768 
00769   // aim for smaller dimension equal to target, larger dimension larger by Ly/Lx or Lx/Ly proportion
00770   const double  yTOx = Ly / Lx;
00771   if (Ly > Lx) {
00772     x = target_size; 
00773     y = (PetscInt) ((double)target_size * yTOx); 
00774   } else {
00775     y = target_size; 
00776     x = (PetscInt) ((double)target_size / yTOx);
00777   }
00778   
00779   // if either dimension is larger than twice the target, shrink appropriately
00780   if (x > 2 * target_size) {
00781     y = (PetscInt) ( (double)(y) * (2.0 * (double)target_size / (double)(x)) );
00782     x = 2 * target_size;
00783   } else if (y > 2 * target_size) {
00784     x = (PetscInt) ( (double)(x) * (2.0 * (double)target_size / (double)(y)) );
00785     y = 2 * target_size;
00786   }
00787   
00788   // make sure minimum dimension is sufficient to see
00789   if (x < 20)   x = 20;
00790   if (y < 20)   y = 20;
00791   return 0;
00792 }
00793 
00795 PetscErrorCode IceModelVec::has_nan() {
00796   PetscErrorCode ierr;
00797   PetscReal tmp;
00798 
00799   ierr = norm(NORM_INFINITY, tmp); CHKERRQ(ierr);
00800 
00801   if ( gsl_isnan(tmp) ) {
00802     //    SETERRQ1(1, "IceModelVec %s has NANs", name.c_str());
00803     PetscPrintf(grid->com, "IceModelVec %s has uninitialized grid points (or NANs)\n", name.c_str());
00804   }
00805 
00806   return 0;
00807 }
00808 
00809 /********* IceModelVec3 and IceModelVec3Bedrock: SEE SEPARATE FILE  iceModelVec3.cc    **********/
00810 
00811 /********* IceModelVec2: SEE SEPARATE FILE  iceModelVec2.cc    **********/