v0.13.0
SnesCtx.cpp
Go to the documentation of this file.
1 /* This file is part of MoFEM.
2  * MoFEM is free software: you can redistribute it and/or modify it under
3  * the terms of the GNU Lesser General Public License as published by the
4  * Free Software Foundation, either version 3 of the License, or (at your
5  * option) any later version.
6  *
7  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10  * License for more details.
11  *
12  * You should have received a copy of the GNU Lesser General Public
13  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14 
15 namespace MoFEM {
16 
21  preProcess_Mat = snes_ctx.preProcess_Mat;
23  preProcess_Rhs = snes_ctx.preProcess_Rhs;
26 }
27 
30  loops_to_do_Mat.clear();
31  loops_to_do_Rhs.clear();
32  preProcess_Mat.clear();
33  postProcess_Mat.clear();
34  preProcess_Rhs.clear();
35  postProcess_Rhs.clear();
37 }
38 
39 PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx) {
40  SnesCtx *snes_ctx = (SnesCtx *)ctx;
41  // PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
43  PetscLogEventBegin(snes_ctx->MOFEM_EVENT_SnesRhs, 0, 0, 0, 0);
44  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
45  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
46  if (snes_ctx->vErify) {
47  // Verify finite elements, check for not a number
48  CHKERR VecAssemblyBegin(f);
49  CHKERR VecAssemblyEnd(f);
50  MPI_Comm comm = PetscObjectComm((PetscObject)f);
51  PetscSynchronizedPrintf(comm, "SNES Verify x\n");
52  const Problem *prb_ptr;
53  CHKERR snes_ctx->mField.get_problem(snes_ctx->problemName, &prb_ptr);
54  CHKERR snes_ctx->mField.getInterface<Tools>()->checkVectorForNotANumber(
55  prb_ptr, COL, x);
56  }
57  CHKERR snes_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
58  snes_ctx->problemName, COL, x, INSERT_VALUES, SCATTER_REVERSE);
59 
60  auto zero_ghost_vec = [](Vec g) {
62  Vec l;
63  CHKERR VecGhostGetLocalForm(g, &l);
64  double *a;
65  CHKERR VecGetArray(l, &a);
66  int s;
67  CHKERR VecGetLocalSize(l, &s);
68  for (int i = 0; i != s; ++i)
69  a[i] = 0;
70  CHKERR VecRestoreArray(l, &a);
71  CHKERR VecGhostRestoreLocalForm(g, &l);
73  };
74  CHKERR zero_ghost_vec(f);
75 
76  snes_ctx->vecAssembleSwitch = boost::movelib::make_unique<bool>(true);
77 
78  auto set = [&](auto &fe) {
79  fe.snes = snes;
80  fe.snes_x = x;
81  fe.snes_f = f;
82  fe.snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
83  fe.ksp_ctx = KspMethod::CTX_SETFUNCTION;
84  fe.data_ctx = PetscData::CtxSetF | PetscData::CtxSetX;
85  };
86 
87  auto unset = [&](auto &fe) {
88  fe.snes_ctx = SnesMethod::CTX_SNESNONE;
89  fe.ksp_ctx = KspMethod::CTX_KSPNONE;
90  fe.data_ctx = PetscData::CtxSetNone;
91  };
92 
93  for (auto &bit : snes_ctx->preProcess_Rhs) {
94  bit->vecAssembleSwitch = boost::move(snes_ctx->vecAssembleSwitch);
95  set(*bit);
97  snes_ctx->problemName, *bit);
98  unset(*bit);
99  snes_ctx->vecAssembleSwitch = boost::move(bit->vecAssembleSwitch);
100  }
101 
102  auto cache_ptr = boost::make_shared<CacheTuple>();
103  CHKERR snes_ctx->mField.cache_problem_entities(snes_ctx->problemName,
104  cache_ptr);
105 
106  for (auto &lit : snes_ctx->loops_to_do_Rhs) {
107  lit.second->vecAssembleSwitch = boost::move(snes_ctx->vecAssembleSwitch);
108  set(*lit.second);
110  snes_ctx->problemName, lit.first, *lit.second, nullptr, snes_ctx->bH,
111  cache_ptr);
112  unset(*lit.second);
113  if (snes_ctx->vErify) {
114  // Verify finite elements, check for not a number
115  CHKERR VecAssemblyBegin(f);
116  CHKERR VecAssemblyEnd(f);
117  MPI_Comm comm = PetscObjectComm((PetscObject)f);
118  PetscSynchronizedPrintf(comm, "SNES Verify f FE < %s >\n",
119  lit.first.c_str());
120  const Problem *prb_ptr;
121  CHKERR snes_ctx->mField.get_problem(snes_ctx->problemName, &prb_ptr);
122  CHKERR snes_ctx->mField.getInterface<Tools>()->checkVectorForNotANumber(
123  prb_ptr, ROW, f);
124  }
125 
126  snes_ctx->vecAssembleSwitch = boost::move(lit.second->vecAssembleSwitch);
127  }
128 
129  for (auto &bit : snes_ctx->postProcess_Rhs) {
130  bit->vecAssembleSwitch = boost::move(snes_ctx->vecAssembleSwitch);
131  set(*bit);
133  snes_ctx->problemName, *bit);
134  unset(*bit);
135  snes_ctx->vecAssembleSwitch = boost::move(bit->vecAssembleSwitch);
136  }
137 
138  if (snes_ctx->vecAssembleSwitch) {
139  CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
140  CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
141  CHKERR VecAssemblyBegin(f);
142  CHKERR VecAssemblyEnd(f);
143  }
144  PetscLogEventEnd(snes_ctx->MOFEM_EVENT_SnesRhs, 0, 0, 0, 0);
146 }
147 
148 PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx) {
149  SnesCtx *snes_ctx = (SnesCtx *)ctx;
150  // PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
152  PetscLogEventBegin(snes_ctx->MOFEM_EVENT_SnesMat, 0, 0, 0, 0);
153  if (snes_ctx->zeroPreCondMatrixB)
154  CHKERR MatZeroEntries(B);
155 
156  snes_ctx->matAssembleSwitch = boost::movelib::make_unique<bool>(true);
157 
158  auto set = [&](auto &fe) {
159  fe.snes = snes;
160  fe.snes_x = x;
161  fe.snes_A = A;
162  fe.snes_B = B;
163  fe.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
164  fe.ksp_ctx = KspMethod::CTX_OPERATORS;
166  };
167 
168  auto unset = [&](auto &fe) {
169  fe.snes_ctx = SnesMethod::CTX_SNESNONE;
170  fe.ksp_ctx = KspMethod::CTX_KSPNONE;
171  fe.data_ctx = PetscData::CtxSetNone;
172  };
173 
174 
175  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
176  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
177  CHKERR snes_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
178  snes_ctx->problemName, COL, x, INSERT_VALUES, SCATTER_REVERSE);
179  for (auto &bit : snes_ctx->preProcess_Mat) {
180  bit->matAssembleSwitch = boost::move(snes_ctx->matAssembleSwitch);
181  set(*bit);
183  snes_ctx->problemName, *bit);
184  unset(*bit);
185  snes_ctx->matAssembleSwitch = boost::move(bit->matAssembleSwitch);
186  }
187 
188  auto cache_ptr = boost::make_shared<CacheTuple>();
189  CHKERR snes_ctx->mField.cache_problem_entities(snes_ctx->problemName,
190  cache_ptr);
191 
192  for (auto &lit : snes_ctx->loops_to_do_Mat) {
193  lit.second->matAssembleSwitch = boost::move(snes_ctx->matAssembleSwitch);
194  set(*lit.second);
196  snes_ctx->problemName, lit.first, *(lit.second), nullptr, snes_ctx->bH,
197  cache_ptr);
198  unset(*lit.second);
199  snes_ctx->matAssembleSwitch = boost::move(lit.second->matAssembleSwitch);
200  }
201 
202  for (auto &bit : snes_ctx->postProcess_Mat) {
203  bit->matAssembleSwitch = boost::move(snes_ctx->matAssembleSwitch);
204  set(*bit);
206  snes_ctx->problemName, *bit);
207  unset(*bit);
208  snes_ctx->matAssembleSwitch = boost::move(bit->matAssembleSwitch);
209  }
210 
211  if (*snes_ctx->matAssembleSwitch) {
212  CHKERR MatAssemblyBegin(B, snes_ctx->typeOfAssembly);
213  CHKERR MatAssemblyEnd(B, snes_ctx->typeOfAssembly);
214  }
215  PetscLogEventEnd(snes_ctx->MOFEM_EVENT_SnesMat, 0, 0, 0, 0);
217 }
218 
219 MoFEMErrorCode SnesMoFEMSetAssemblyType(SNES snes, MatAssemblyType type) {
220  SnesCtx *snes_ctx;
221  // PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
223  CHKERR SNESGetApplicationContext(snes, &snes_ctx);
224  snes_ctx->typeOfAssembly = type;
226 }
227 
229  SnesCtx *snes_ctx;
231  CHKERR SNESGetApplicationContext(snes, &snes_ctx);
232  snes_ctx->bH = bh;
234 }
235 
236 } // namespace MoFEM
constexpr double a
@ COL
Definition: definitions.h:136
@ ROW
Definition: definitions.h:136
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
Definition: definitions.h:110
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
virtual MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition: SnesCtx.cpp:148
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition: SnesCtx.cpp:39
MoFEMErrorCode SnesMoFEMSetAssemblyType(SNES snes, MatAssemblyType type)
Set assembly type at the end of SnesMat.
Definition: SnesCtx.cpp:219
MoFEMErrorCode SnesMoFEMSetBehavior(SNES snes, MoFEMTypes bh)
Set behavior if finite element in sequence does not exist.
Definition: SnesCtx.cpp:228
double A
constexpr double g
virtual MoFEMErrorCode cache_problem_entities(const std::string prb_name, CacheTupleWeakPtr cache_ptr)=0
Cache variables.
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
static constexpr Switches CtxSetA
Definition: LoopMethods.hpp:47
static constexpr Switches CtxSetX
Definition: LoopMethods.hpp:49
static constexpr Switches CtxSetNone
Definition: LoopMethods.hpp:45
static constexpr Switches CtxSetF
Definition: LoopMethods.hpp:46
static constexpr Switches CtxSetB
Definition: LoopMethods.hpp:48
keeps basic data about problem
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:27
PetscLogEvent MOFEM_EVENT_SnesMat
Log events to assemble tangent matrix.
Definition: SnesCtx.hpp:137
BasicMethodsSequence preProcess_Rhs
Sequence of methods run before residual is assembled.
Definition: SnesCtx.hpp:53
BasicMethodsSequence postProcess_Mat
Definition: SnesCtx.hpp:50
MoFEMTypes bH
If set to MF_EXIST check if element exist, default MF_EXIST.
Definition: SnesCtx.hpp:34
BasicMethodsSequence postProcess_Rhs
Sequence of methods run after residual is assembled.
Definition: SnesCtx.hpp:55
MoFEMErrorCode copyLoops(const SnesCtx &snes_ctx)
Copy sequences from other SNES contex.
Definition: SnesCtx.cpp:17
MatAssemblyType typeOfAssembly
type of assembly at the end
Definition: SnesCtx.hpp:37
bool zeroPreCondMatrixB
Definition: SnesCtx.hpp:35
MoFEMErrorCode clearLoops()
Clear loops.
Definition: SnesCtx.cpp:28
BasicMethodsSequence preProcess_Mat
Definition: SnesCtx.hpp:48
bool vErify
If true verify vector.
Definition: SnesCtx.hpp:38
FEMethodsSequence loops_to_do_Mat
Definition: SnesCtx.hpp:44
std::string problemName
problem name
Definition: SnesCtx.hpp:32
FEMethodsSequence loops_to_do_Rhs
Definition: SnesCtx.hpp:46
boost::movelib::unique_ptr< bool > matAssembleSwitch
Definition: SnesCtx.hpp:135
MoFEM::Interface & mField
database Interface
Definition: SnesCtx.hpp:29
boost::movelib::unique_ptr< bool > vecAssembleSwitch
Definition: SnesCtx.hpp:134
PetscLogEvent MOFEM_EVENT_SnesRhs
Log events to assemble residual.
Definition: SnesCtx.hpp:136
Auxiliary tools.
Definition: Tools.hpp:28
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:33