PISM, A Parallel Ice Sheet Model  stable v2.1.1 committed by Constantine Khrulev on 2024-12-04 13:36:58 -0900
SNESProblem.hh
Go to the documentation of this file.
1 // Copyright (C) 2011, 2012, 2014, 2015, 2016, 2017, 2022, 2024 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 
20 #ifndef PISM_SNESPROBLEM_H
21 #define PISM_SNESPROBLEM_H
22 
23 #include "pism/util/Grid.hh" // inline implementation in the header uses Grid
24 #include "pism/util/Vector2d.hh" // to get Vector2
25 #include "pism/util/petscwrappers/SNES.hh"
26 #include "pism/util/Logger.hh"
27 
28 namespace pism {
29 
30 template<int DOF, class U> class SNESProblem {
31 public:
32  SNESProblem(std::shared_ptr<const Grid> g);
33 
34  virtual ~SNESProblem();
35 
36  virtual void solve();
37 
38  virtual const std::string& name();
39 
40  virtual Vec solution()
41  {
42  return m_X;
43  }
44 
45 protected:
46 
47  virtual void compute_local_function(DMDALocalInfo *info, const U **xg, U **yg) = 0;
48  virtual void compute_local_jacobian(DMDALocalInfo *info, const U **x, Mat B) = 0;
49 
50  std::shared_ptr<const Grid> m_grid;
51 
55 
56 private:
57 
58  struct CallbackData {
59  DM da;
61  };
62 
64 
65  static PetscErrorCode function_callback(DMDALocalInfo *info, const U **x, U **f,
66  CallbackData *);
67  static PetscErrorCode jacobian_callback(DMDALocalInfo *info, const U **x, Mat B,
68  CallbackData *);
69 };
70 
73 
74 template<int DOF, class U>
75 PetscErrorCode SNESProblem<DOF,U>::function_callback(DMDALocalInfo *info,
76  const U **x, U **f,
78  try {
79  cb->solver->compute_local_function(info,x,f);
80  } catch (...) {
81  MPI_Comm com = MPI_COMM_SELF;
82  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)cb->da, &com); CHKERRQ(ierr);
84  SETERRQ(com, 1, "A PISM callback failed");
85  }
86  return 0;
87 }
88 
89 template<int DOF, class U>
90 PetscErrorCode SNESProblem<DOF,U>::jacobian_callback(DMDALocalInfo *info,
91  const U **x, Mat J,
93  try {
94  cb->solver->compute_local_jacobian(info, x, J);
95  } catch (...) {
96  MPI_Comm com = MPI_COMM_SELF;
97  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)cb->da, &com); CHKERRQ(ierr);
99  SETERRQ(com, 1, "A PISM callback failed");
100  }
101  return 0;
102 }
103 
104 template<int DOF, class U>
105 SNESProblem<DOF, U>::SNESProblem(std::shared_ptr<const Grid> g)
106  : m_grid(g) {
107 
108  PetscErrorCode ierr;
109 
110  int stencil_width=1;
111  m_DA = m_grid->get_dm(DOF, stencil_width);
112 
113  ierr = DMCreateGlobalVector(*m_DA, m_X.rawptr());
114  PISM_CHK(ierr, "DMCreateGlobalVector");
115 
116  ierr = SNESCreate(m_grid->com, m_snes.rawptr());
117  PISM_CHK(ierr, "SNESCreate");
118 
119  // Set the SNES callbacks to call into our compute_local_function and compute_local_jacobian
121  m_callbackData.solver = this;
122 
123  ierr = DMDASNESSetFunctionLocal(*m_DA, INSERT_VALUES,
124 #if PETSC_VERSION_LT(3,21,0)
125  (DMDASNESFunction)SNESProblem<DOF, U>::function_callback,
126 #else
127  (DMDASNESFunctionFn*)SNESProblem<DOF, U>::function_callback,
128 #endif
129  &m_callbackData);
130  PISM_CHK(ierr, "DMDASNESSetFunctionLocal");
131 
132  ierr = DMDASNESSetJacobianLocal(*m_DA,
133 #if PETSC_VERSION_LT(3,21,0)
134  (DMDASNESJacobian)SNESProblem<DOF, U>::jacobian_callback,
135 #else
136  (DMDASNESJacobianFn*)SNESProblem<DOF, U>::jacobian_callback,
137 #endif
138  &m_callbackData);
139  PISM_CHK(ierr, "DMDASNESSetJacobianLocal");
140 
141  ierr = DMSetMatType(*m_DA, "baij");
142  PISM_CHK(ierr, "DMSetMatType");
143 
144  ierr = DMSetApplicationContext(*m_DA, &m_callbackData);
145  PISM_CHK(ierr, "DMSetApplicationContext");
146 
147  ierr = SNESSetDM(m_snes, *m_DA);
148  PISM_CHK(ierr, "SNESSetDM");
149 
150  ierr = SNESSetFromOptions(m_snes);
151  PISM_CHK(ierr, "SNESSetFromOptions");
152 }
153 
154 template<int DOF, class U>
156  // empty
157 }
158 
159 template<int DOF, class U>
160 const std::string& SNESProblem<DOF,U>::name() {
161  return "UnnamedProblem";
162 }
163 
164 template<int DOF, class U>
166  PetscErrorCode ierr;
167 
168  // Solve:
169  ierr = SNESSolve(m_snes,NULL,m_X); PISM_CHK(ierr, "SNESSolve");
170 
171  // See if it worked.
172  SNESConvergedReason reason;
173  ierr = SNESGetConvergedReason(m_snes, &reason); PISM_CHK(ierr, "SNESGetConvergedReason");
174  if (reason < 0) {
175  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "SNESProblem %s solve failed to converge (SNES reason %s)",
176  name().c_str(), SNESConvergedReasons[reason]);
177  }
178 
179  m_grid->ctx()->log()->message(1, "SNESProblem %s converged (SNES reason %s)\n",
180  name().c_str(), SNESConvergedReasons[reason]);
181 }
182 
183 } // end of namespace pism
184 
185 #endif // PISM_SNESPROBLEM_H
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
SNESProblem(std::shared_ptr< const Grid > g)
Definition: SNESProblem.hh:105
CallbackData m_callbackData
Definition: SNESProblem.hh:63
virtual const std::string & name()
Definition: SNESProblem.hh:160
petsc::DM::Ptr m_DA
Definition: SNESProblem.hh:54
virtual void compute_local_function(DMDALocalInfo *info, const U **xg, U **yg)=0
static PetscErrorCode jacobian_callback(DMDALocalInfo *info, const U **x, Mat B, CallbackData *)
Definition: SNESProblem.hh:90
petsc::Vec m_X
Definition: SNESProblem.hh:52
virtual void compute_local_jacobian(DMDALocalInfo *info, const U **x, Mat B)=0
virtual Vec solution()
Definition: SNESProblem.hh:40
std::shared_ptr< const Grid > m_grid
Definition: SNESProblem.hh:50
virtual void solve()
Definition: SNESProblem.hh:165
petsc::SNES m_snes
Definition: SNESProblem.hh:53
virtual ~SNESProblem()
Definition: SNESProblem.hh:155
static PetscErrorCode function_callback(DMDALocalInfo *info, const U **x, U **f, CallbackData *)
Definition: SNESProblem.hh:75
T * rawptr()
Definition: Wrapper.hh:39
std::shared_ptr< Wrapper > Ptr
Definition: Wrapper.hh:30
#define PISM_CHK(errcode, name)
#define PISM_ERROR_LOCATION
static const double g
Definition: exactTestP.cc:36
SNESProblem< 1, double > SNESScalarProblem
Definition: SNESProblem.hh:71
void handle_fatal_errors(MPI_Comm com)
SNESProblem< 2, Vector2d > SNESVectorProblem
Definition: SNESProblem.hh:72
const int J[]
Definition: ssafd_code.cc:34
SNESProblem< DOF, U > * solver
Definition: SNESProblem.hh:60