v0.14.0
SnesCtx.cpp
Go to the documentation of this file.
1 
2 
3 namespace MoFEM {
4 
11  preProcess_Rhs = snes_ctx.preProcess_Rhs;
14 }
15 
18  loops_to_do_Mat.clear();
19  loops_to_do_Rhs.clear();
20  preProcess_Mat.clear();
21  postProcess_Mat.clear();
22  preProcess_Rhs.clear();
23  postProcess_Rhs.clear();
25 }
26 
27 PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx) {
28  SnesCtx *snes_ctx = (SnesCtx *)ctx;
29  // PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
31  PetscLogEventBegin(snes_ctx->MOFEM_EVENT_SnesRhs, 0, 0, 0, 0);
32  Vec dx;
33  CHKERR SNESGetSolutionUpdate(snes, &dx);
34  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
35  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
36  if (snes_ctx->vErify) {
37  // Verify finite elements, check for not a number
38  CHKERR VecAssemblyBegin(f);
39  CHKERR VecAssemblyEnd(f);
40  MPI_Comm comm = PetscObjectComm((PetscObject)f);
41  PetscSynchronizedPrintf(comm, "SNES Verify x\n");
42  const Problem *prb_ptr;
43  CHKERR snes_ctx->mField.get_problem(snes_ctx->problemName, &prb_ptr);
44  CHKERR snes_ctx->mField.getInterface<Tools>()->checkVectorForNotANumber(
45  prb_ptr, COL, x);
46  }
47  CHKERR snes_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
48  snes_ctx->problemName, COL, x, INSERT_VALUES, SCATTER_REVERSE);
49 
50  auto zero_ghost_vec = [](Vec g) {
52  Vec l;
53  CHKERR VecGhostGetLocalForm(g, &l);
54  double *a;
55  CHKERR VecGetArray(l, &a);
56  int s;
57  CHKERR VecGetLocalSize(l, &s);
58  for (int i = 0; i != s; ++i)
59  a[i] = 0;
60  CHKERR VecRestoreArray(l, &a);
61  CHKERR VecGhostRestoreLocalForm(g, &l);
63  };
64  CHKERR zero_ghost_vec(f);
65 
66  snes_ctx->vecAssembleSwitch = boost::movelib::make_unique<bool>(true);
67  auto cache_ptr = boost::make_shared<CacheTuple>();
68  CHKERR snes_ctx->mField.cache_problem_entities(snes_ctx->problemName,
69  cache_ptr);
70 
71  auto set = [&](auto &fe) {
72  fe.snes = snes;
73  fe.snes_x = x;
74  fe.snes_dx = dx;
75  fe.snes_f = f;
76  fe.snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
77  fe.ksp_ctx = KspMethod::CTX_SETFUNCTION;
79 
80  CHKERR SNESGetKSP(snes, &fe.ksp);
81 
82  fe.cacheWeakPtr = cache_ptr;
83  };
84 
85  auto unset = [&](auto &fe) {
86  fe.snes_ctx = SnesMethod::CTX_SNESNONE;
87  fe.ksp_ctx = KspMethod::CTX_KSPNONE;
88  fe.data_ctx = PetscData::CtxSetNone;
89  };
90 
91  for (auto &bit : snes_ctx->preProcess_Rhs) {
92  bit->vecAssembleSwitch = boost::move(snes_ctx->vecAssembleSwitch);
93  set(*bit);
95  snes_ctx->problemName, *bit);
96  unset(*bit);
97  snes_ctx->vecAssembleSwitch = boost::move(bit->vecAssembleSwitch);
98  }
99 
100  for (auto &lit : snes_ctx->loops_to_do_Rhs) {
101  lit.second->vecAssembleSwitch = boost::move(snes_ctx->vecAssembleSwitch);
102  set(*lit.second);
104  snes_ctx->problemName, lit.first, *lit.second, nullptr, snes_ctx->bH,
105  cache_ptr);
106  unset(*lit.second);
107  if (snes_ctx->vErify) {
108  // Verify finite elements, check for not a number
109  CHKERR VecAssemblyBegin(f);
110  CHKERR VecAssemblyEnd(f);
111  MPI_Comm comm = PetscObjectComm((PetscObject)f);
112  PetscSynchronizedPrintf(comm, "SNES Verify f FE < %s >\n",
113  lit.first.c_str());
114  const Problem *prb_ptr;
115  CHKERR snes_ctx->mField.get_problem(snes_ctx->problemName, &prb_ptr);
116  CHKERR snes_ctx->mField.getInterface<Tools>()->checkVectorForNotANumber(
117  prb_ptr, ROW, f);
118  }
119 
120  snes_ctx->vecAssembleSwitch = boost::move(lit.second->vecAssembleSwitch);
121  }
122 
123  for (auto &bit : snes_ctx->postProcess_Rhs) {
124  bit->vecAssembleSwitch = boost::move(snes_ctx->vecAssembleSwitch);
125  set(*bit);
127  snes_ctx->problemName, *bit);
128  unset(*bit);
129  snes_ctx->vecAssembleSwitch = boost::move(bit->vecAssembleSwitch);
130  }
131 
132  if (*(snes_ctx->vecAssembleSwitch)) {
133  CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
134  CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
135  CHKERR VecAssemblyBegin(f);
136  CHKERR VecAssemblyEnd(f);
137  }
138  PetscLogEventEnd(snes_ctx->MOFEM_EVENT_SnesRhs, 0, 0, 0, 0);
140 }
141 
142 PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx) {
143  SnesCtx *snes_ctx = (SnesCtx *)ctx;
144  // PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
146  PetscLogEventBegin(snes_ctx->MOFEM_EVENT_SnesMat, 0, 0, 0, 0);
147  if (snes_ctx->zeroPreCondMatrixB)
148  CHKERR MatZeroEntries(B);
149 
150  snes_ctx->matAssembleSwitch = boost::movelib::make_unique<bool>(true);
151  auto cache_ptr = boost::make_shared<CacheTuple>();
152  CHKERR snes_ctx->mField.cache_problem_entities(snes_ctx->problemName,
153  cache_ptr);
154  Vec dx;
155  CHKERR SNESGetSolutionUpdate(snes, &dx);
156 
157  auto set = [&](auto &fe) {
158  fe.snes = snes;
159  fe.snes_x = x;
160  fe.snes_dx = dx;
161  fe.snes_A = A;
162  fe.snes_B = B;
163  fe.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
164  fe.ksp_ctx = KspMethod::CTX_OPERATORS;
167 
168  CHKERR SNESGetKSP(snes, &fe.ksp);
169 
170  fe.cacheWeakPtr = cache_ptr;
171  };
172 
173  auto unset = [&](auto &fe) {
174  fe.snes_ctx = SnesMethod::CTX_SNESNONE;
175  fe.ksp_ctx = KspMethod::CTX_KSPNONE;
176  fe.data_ctx = PetscData::CtxSetNone;
177  };
178 
179 
180  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
181  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
182  CHKERR snes_ctx->mField.getInterface<VecManager>()->setLocalGhostVector(
183  snes_ctx->problemName, COL, x, INSERT_VALUES, SCATTER_REVERSE);
184  for (auto &bit : snes_ctx->preProcess_Mat) {
185  bit->matAssembleSwitch = boost::move(snes_ctx->matAssembleSwitch);
186  set(*bit);
188  snes_ctx->problemName, *bit);
189  unset(*bit);
190  snes_ctx->matAssembleSwitch = boost::move(bit->matAssembleSwitch);
191  }
192 
193 
194 
195  for (auto &lit : snes_ctx->loops_to_do_Mat) {
196  lit.second->matAssembleSwitch = boost::move(snes_ctx->matAssembleSwitch);
197  set(*lit.second);
199  snes_ctx->problemName, lit.first, *(lit.second), nullptr, snes_ctx->bH,
200  cache_ptr);
201  unset(*lit.second);
202  snes_ctx->matAssembleSwitch = boost::move(lit.second->matAssembleSwitch);
203  }
204 
205  for (auto &bit : snes_ctx->postProcess_Mat) {
206  bit->matAssembleSwitch = boost::move(snes_ctx->matAssembleSwitch);
207  set(*bit);
209  snes_ctx->problemName, *bit);
210  unset(*bit);
211  snes_ctx->matAssembleSwitch = boost::move(bit->matAssembleSwitch);
212  }
213 
214  if (*snes_ctx->matAssembleSwitch) {
215  CHKERR MatAssemblyBegin(B, snes_ctx->typeOfAssembly);
216  CHKERR MatAssemblyEnd(B, snes_ctx->typeOfAssembly);
217  }
218  PetscLogEventEnd(snes_ctx->MOFEM_EVENT_SnesMat, 0, 0, 0, 0);
220 }
221 
222 MoFEMErrorCode SnesMoFEMSetAssemblyType(SNES snes, MatAssemblyType type) {
223  SnesCtx *snes_ctx;
224  // PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
226  CHKERR SNESGetApplicationContext(snes, &snes_ctx);
227  snes_ctx->typeOfAssembly = type;
229 }
230 
232  SnesCtx *snes_ctx;
234  CHKERR SNESGetApplicationContext(snes, &snes_ctx);
235  snes_ctx->bH = bh;
237 }
238 
239 MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm,
240  SnesCtx *snes_ctx) {
242  auto &m_field = snes_ctx->mField;
243  auto problem_ptr = m_field.get_problem(snes_ctx->problemName);
244  auto fields_ptr = m_field.get_fields();
245  auto dofs = problem_ptr->numeredRowDofsPtr;
246 
247  std::vector<double> lnorms(fields_ptr->size(), 0),
248  norms(fields_ptr->size(), 0);
249 
250  Vec res;
251  CHKERR SNESGetFunction(snes, &res, NULL, NULL);
252 
253  const double *r;
254  CHKERR VecGetArrayRead(res, &r);
255  {
256  int f = 0;
257  for (auto fi : *fields_ptr) {
258  const auto lo_uid = FieldEntity::getLoBitNumberUId(fi->bitNumber);
259  const auto hi_uid = FieldEntity::getHiBitNumberUId(fi->bitNumber);
260  const auto hi = dofs->get<Unique_mi_tag>().upper_bound(hi_uid);
261  for (auto lo = dofs->get<Unique_mi_tag>().lower_bound(lo_uid); lo != hi;
262  ++lo) {
263  const DofIdx loc_idx = (*lo)->getPetscLocalDofIdx();
264  if (loc_idx >= 0 && loc_idx < problem_ptr->nbLocDofsRow) {
265  lnorms[f] += PetscRealPart(PetscSqr(r[loc_idx]));
266  }
267  }
268  ++f;
269  }
270  }
271  CHKERR VecRestoreArrayRead(res, &r);
272 
273  MPIU_Allreduce(&*lnorms.begin(), &*norms.begin(), lnorms.size(), MPIU_REAL,
274  MPIU_SUM, PetscObjectComm((PetscObject)snes));
275 
276  std::stringstream s;
277  int tl;
278  CHKERR PetscObjectGetTabLevel((PetscObject)snes, &tl);
279  for (auto t = 0; t != tl; ++t)
280  s << " ";
281  s << its << " Function norm " << boost::format("%14.12e") % (double)fgnorm
282  << " [";
283  {
284  int f = 0;
285  for (auto fi : *fields_ptr) {
286  if (f > 0)
287  s << ", ";
288  s << boost::format("%14.12e") % (double)PetscSqrtReal(norms[f]);
289  ++f;
290  }
291  s << "]";
292  }
293 
294  MOFEM_LOG("SNES_WORLD", Sev::inform) << s.str();
295 
297 }
298 
299 } // namespace MoFEM
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::KspMethod::CTX_KSPNONE
@ CTX_KSPNONE
Definition: LoopMethods.hpp:76
g
constexpr double g
Definition: shallow_wave.cpp:63
MoFEM::CoreInterface::problem_basic_method_postProcess
virtual MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
MoFEM::CoreInterface::loop_finite_elements
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.
MoFEM::SnesCtx::MOFEM_EVENT_SnesMat
PetscLogEvent MOFEM_EVENT_SnesMat
Log events to assemble tangent matrix.
Definition: SnesCtx.hpp:158
MoFEM::SnesCtx::vecAssembleSwitch
boost::movelib::unique_ptr< bool > vecAssembleSwitch
Definition: SnesCtx.hpp:155
MoFEM::PetscData::CtxSetB
static constexpr Switches CtxSetB
Definition: LoopMethods.hpp:39
MoFEM::Tools
Auxiliary tools.
Definition: Tools.hpp:19
MoFEM::PetscData::CtxSetA
static constexpr Switches CtxSetA
Definition: LoopMethods.hpp:38
MoFEM::SnesCtx::vErify
bool vErify
If true verify vector.
Definition: SnesCtx.hpp:24
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::CoreInterface::get_problem
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::SnesCtx::postProcess_Mat
BasicMethodsSequence postProcess_Mat
Definition: SnesCtx.hpp:36
MoFEM::KspMethod::CTX_SETFUNCTION
@ CTX_SETFUNCTION
Definition: LoopMethods.hpp:76
MoFEM::SnesMethod::CTX_SNESSETJACOBIAN
@ CTX_SNESSETJACOBIAN
Definition: LoopMethods.hpp:110
sdf.r
int r
Definition: sdf.py:8
MoFEM::SnesMethod::CTX_SNESSETFUNCTION
@ CTX_SNESSETFUNCTION
Definition: LoopMethods.hpp:110
ROW
@ ROW
Definition: definitions.h:136
MoFEM::SnesCtx::loops_to_do_Rhs
FEMethodsSequence loops_to_do_Rhs
Definition: SnesCtx.hpp:32
MoFEM::CoreInterface::problem_basic_method_preProcess
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::SnesCtx::copyLoops
MoFEMErrorCode copyLoops(const SnesCtx &snes_ctx)
Copy sequences from other SNES contex.
Definition: SnesCtx.cpp:5
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::SnesCtx::zeroPreCondMatrixB
bool zeroPreCondMatrixB
Definition: SnesCtx.hpp:21
MoFEM::SnesCtx::mField
MoFEM::Interface & mField
database Interface
Definition: SnesCtx.hpp:15
MoFEM::SnesCtx::clearLoops
MoFEMErrorCode clearLoops()
Clear loops.
Definition: SnesCtx.cpp:16
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::SnesCtx::typeOfAssembly
MatAssemblyType typeOfAssembly
type of assembly at the end
Definition: SnesCtx.hpp:23
MoFEM::SnesMethod::CTX_SNESNONE
@ CTX_SNESNONE
Definition: LoopMethods.hpp:110
double
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
convert.type
type
Definition: convert.py:64
MoFEM::SnesCtx::problemName
std::string problemName
problem name
Definition: SnesCtx.hpp:18
MoFEM::FieldEntity::getLoBitNumberUId
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
Definition: FieldEntsMultiIndices.hpp:222
MoFEM::SnesCtx::postProcess_Rhs
BasicMethodsSequence postProcess_Rhs
Sequence of methods run after residual is assembled.
Definition: SnesCtx.hpp:41
COL
@ COL
Definition: definitions.h:136
MoFEM::SnesCtx::MOFEM_EVENT_SnesRhs
PetscLogEvent MOFEM_EVENT_SnesRhs
Log events to assemble residual.
Definition: SnesCtx.hpp:157
MoFEM::CoreInterface::cache_problem_entities
virtual MoFEMErrorCode cache_problem_entities(const std::string prb_name, CacheTupleWeakPtr cache_ptr)=0
Cache variables.
MoFEM::SnesMoFEMSetBehavior
MoFEMErrorCode SnesMoFEMSetBehavior(SNES snes, MoFEMTypes bh)
Set behavior if finite element in sequence does not exist.
Definition: SnesCtx.cpp:231
MoFEM::SnesCtx::preProcess_Rhs
BasicMethodsSequence preProcess_Rhs
Sequence of methods run before residual is assembled.
Definition: SnesCtx.hpp:39
MoFEM::PetscData::CtxSetNone
static constexpr Switches CtxSetNone
Definition: LoopMethods.hpp:36
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
MoFEM::Problem::numeredRowDofsPtr
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
Definition: ProblemsMultiIndices.hpp:73
MoFEM::PetscData::CtxSetX
static constexpr Switches CtxSetX
Definition: LoopMethods.hpp:40
MoFEM::SnesCtx::loops_to_do_Mat
FEMethodsSequence loops_to_do_Mat
Definition: SnesCtx.hpp:30
MoFEM::SnesRhs
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:27
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::SnesCtx::matAssembleSwitch
boost::movelib::unique_ptr< bool > matAssembleSwitch
Definition: SnesCtx.hpp:156
MoFEM::PetscData::CtxSetDX
static constexpr Switches CtxSetDX
Definition: LoopMethods.hpp:41
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
MoFEM::SnesCtx::preProcess_Mat
BasicMethodsSequence preProcess_Mat
Definition: SnesCtx.hpp:34
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:64
MoFEM::SnesCtx::bH
MoFEMTypes bH
If set to MF_EXIST check if element exist, default MF_EXIST.
Definition: SnesCtx.hpp:20
MoFEM::MoFEMSNESMonitorFields
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:239
MoFEMTypes
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
Definition: definitions.h:110
MoFEM::FieldEntity::getHiBitNumberUId
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
Definition: FieldEntsMultiIndices.hpp:228
MoFEM::KspMethod::CTX_OPERATORS
@ CTX_OPERATORS
Definition: LoopMethods.hpp:76
MoFEM::PetscData::CtxSetF
static constexpr Switches CtxSetF
Definition: LoopMethods.hpp:37
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEM::SnesCtx
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:13
MoFEM::Unique_mi_tag
Definition: TagMultiIndices.hpp:18
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
MoFEM::SnesMat
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:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MoFEM::SnesMoFEMSetAssemblyType
MoFEMErrorCode SnesMoFEMSetAssemblyType(SNES snes, MatAssemblyType type)
Set assembly type at the end of SnesMat.
Definition: SnesCtx.cpp:222
MoFEM::Types::DofIdx
int DofIdx
Index of DOF.
Definition: Types.hpp:18