PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
DataAccess.hh
Go to the documentation of this file.
1 /* Copyright (C) 2020 PISM Authors
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 
20 #ifndef PISM_DATAACCESS_H
21 #define PISM_DATAACCESS_H
22 
23 #include <cassert>
24 #include <petscdmda.h>
25 
26 #include "pism/util/error_handling.hh"
27 
28 namespace pism {
29 
31 /*!
32  * This template class manages access to 2D and 3D Vecs stored in a DM using
33  * `PetscObjectCompose`. Performs cleanup at the end of scope.
34  *
35  * @param[in] da SNES DM for the solution containing 2D and 3D DMs and Vecs
36  * @param[in] dim number of dimensions (2 or 3)
37  * @param[in] type NOT_GHOSTED -- for setting parameters; GHOSTED -- for accessing ghosts
38  * during residual and Jacobian evaluation
39  */
40 template<typename T>
41 class DataAccess {
42 public:
43  DataAccess(DM da, int dim, AccessType type)
44  : m_local(type == GHOSTED) {
45  int ierr;
46 
47  assert(dim == 2 or dim == 3);
48 
49  if (dim == 2) {
50  ierr = setup(da, "2D_DM", "2D_DM_data");
51  } else {
52  ierr = setup(da, "3D_DM", "3D_DM_data");
53  }
54 
55  if (ierr != 0) {
57  "Failed to create an DataAccess instance");
58  }
59  }
60 
61  int setup(DM da, const char *dm_name, const char *vec_name) {
62  PetscErrorCode ierr;
63 
64  m_com = MPI_COMM_SELF;
65  ierr = PetscObjectGetComm((PetscObject)da, &m_com); CHKERRQ(ierr);
66 
67  ierr = PetscObjectQuery((PetscObject)da, dm_name, (PetscObject*)&m_da); CHKERRQ(ierr);
68 
69  if (!m_da) {
70  SETERRQ(m_com, 1, "Failed to get the inner DM");
71  }
72 
73  Vec X;
74  ierr = PetscObjectQuery((PetscObject)da, vec_name, (PetscObject*)&X); CHKERRQ(ierr);
75 
76  if (!X) {
77  SETERRQ(m_com, 1, "Failed to get the inner Vec");
78  }
79 
80  if (m_local) {
81  ierr = DMGetLocalVector(m_da, &m_x); CHKERRQ(ierr);
82 
83  ierr = DMGlobalToLocalBegin(m_da, X, INSERT_VALUES, m_x); CHKERRQ(ierr);
84 
85  ierr = DMGlobalToLocalEnd(m_da, X, INSERT_VALUES, m_x); CHKERRQ(ierr);
86  } else {
87  m_x = X;
88  }
89 
90  ierr = DMDAVecGetArray(m_da, m_x, &m_a); CHKERRQ(ierr);
91 
92  return 0;
93  }
94 
96  try {
97  PetscErrorCode ierr = DMDAVecRestoreArray(m_da, m_x, &m_a);
98  PISM_CHK(ierr, "DMDAVecRestoreArray");
99 
100  if (m_local) {
101  ierr = DMRestoreLocalVector(m_da, &m_x);
102  PISM_CHK(ierr, "DMRestoreLocalVector");
103  }
104  } catch (...) {
106  }
107  }
108 
109  operator T() {
110  return m_a;
111  }
112 private:
113  MPI_Comm m_com;
114  bool m_local;
115  DM m_da;
116  Vec m_x;
117  T m_a;
118 };
119 
120 } // end of namespace pism
121 
122 #endif /* PISM_DATAACCESS_H */
int setup(DM da, const char *dm_name, const char *vec_name)
Definition: DataAccess.hh:61
DataAccess(DM da, int dim, AccessType type)
Definition: DataAccess.hh:43
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
#define PISM_CHK(errcode, name)
#define PISM_ERROR_LOCATION
AccessType
Definition: DataAccess.hh:30
@ GHOSTED
Definition: DataAccess.hh:30
@ NOT_GHOSTED
Definition: DataAccess.hh:30
void handle_fatal_errors(MPI_Comm com)