26 std::vector<GlobalVector*>&
x;
27 std::vector<GlobalVector*>
const&
x_prev;
36 auto context =
static_cast<PetscContext*
>(petsc_context);
38 DBUG(
"PETScNonlinearSolver: residual callback called.");
40 VecCopy(x, context->x[context->process_id]->getRawVector());
44 time_assembly.
start();
45 context->system->assemble(context->x, context->x_prev, context->process_id);
47 INFO(
"[time] Assembly took {} s.", time_assembly.
elapsed());
48 context->system->getResidual(*context->x[context->process_id],
49 *context->x_prev[context->process_id],
51 context->r->finalizeAssembly();
52 context->J->finalizeAssembly();
54 context->system->getJacobian(*context->J);
55 context->system->applyKnownSolutionsNewton(
56 *context->J, *context->r, *context->x[context->process_id]);
58 VecCopy(context->r->getRawVector(), petsc_r);
64 Mat ,
void* petsc_context)
66 DBUG(
"PETScNonlinearSolver: jacobian callback called.");
70 auto context =
static_cast<PetscContext*
>(petsc_context);
71 MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
72 MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
73 MatCopy(context->J->getRawMatrix(), J, DIFFERENT_NONZERO_PATTERN);
87 prefix = prefix +
"_";
95 PetscReal atol, rtol, stol;
97 SNESGetTolerances(
_snes_solver, &atol, &rtol, &stol,
nullptr, &maxf);
98 SNESSetTolerances(
_snes_solver, atol, rtol, stol, maxiter, maxf);
102 PetscOptionsView(
nullptr, PETSC_VIEWER_STDOUT_WORLD);
120 std::vector<GlobalVector*>
const& ,
121 std::vector<GlobalVector*>
const& ,
int const )
129 "Non-equilibrium initial residuum is not implemented for the "
130 "PETScNonlinearSolver.");
134 std::vector<GlobalVector*>& x,
135 std::vector<GlobalVector*>
const& x_prev,
138 std::vector<GlobalVector*>
const&)>
const& ,
139 int const process_id)
141 DBUG(
"PETScNonlinearSolver: solve()");
142 using TimeDiscretizedSystem =
148 DBUG(
"PETScNonlinearSolver: create vectors");
151 system->getMatrixSpecifications(process_id),
_jacobian_id);
153 system->getMatrixSpecifications(process_id),
_residual_id);
160 double time_dirichlet = 0.0;
162 timer_dirichlet.
start();
163 system->computeKnownSolutions(*x[process_id], process_id);
164 system->applyKnownSolutions(*x[process_id]);
165 time_dirichlet += timer_dirichlet.
elapsed();
166 INFO(
"[time] Applying Dirichlet BCs took {} s.", time_dirichlet);
171 system->getMatrixSpecifications(process_id),
_petsc_x_id);
176 DBUG(
"PETScNonlinearSolver: set function");
180 DBUG(
"PETScNonlinearSolver: set jacobian");
182 SNESSetJacobian(
_snes_solver, J_snes.getRawMatrix(), J_snes.getRawMatrix(),
185 std::unique_ptr<GlobalVector> xl =
nullptr;
186 std::unique_ptr<GlobalVector> xu =
nullptr;
190 if ((std::strcmp(snes_type, SNESVINEWTONRSLS) == 0) ||
191 (std::strcmp(snes_type, SNESVINEWTONSSLS) == 0))
194 DBUG(
"PETScNonlinearSolver: set constraints");
196 system->getMatrixSpecifications(process_id));
198 system->getMatrixSpecifications(process_id));
200 system->updateConstraints(*xl, *xu, process_id);
204 SNESVISetVariableBounds(
_snes_solver, xl->getRawVector(),
208 DBUG(
"PETScNonlinearSolver: call SNESSolve");
209 SNESSolve(
_snes_solver,
nullptr, x_snes.getRawVector());
211 SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
213 INFO(
"PETSsSNES convergence reason {}.", reason);
217 INFO(
"PETScSNES used {} iterations.", iterations);
228 return {reason >= 0, iterations};
void INFO(char const *fmt, Args const &... args)
void DBUG(char const *fmt, Args const &... args)
Definition of the RunTime class.
double elapsed() const
Get the elapsed time in seconds.
void start()
Start the timer.
Global vector based on Eigen vector.
virtual GlobalMatrix & getMatrix(std::size_t &id)=0
Get an uninitialized matrix with the given id.
virtual void releaseMatrix(GlobalMatrix const &A)=0
std::size_t _petsc_residual_id
ConvergenceCriterion * _convergence_criterion
NonlinearSolverStatus solve(std::vector< GlobalVector * > &x, std::vector< GlobalVector * > const &x_prev, std::function< void(int, std::vector< GlobalVector * > const &)> const &postIterationCallback, int const process_id) override
std::size_t _petsc_jacobian_id
System * _equation_system
void setEquationSystem(System &eq, ConvergenceCriterion &conv_crit)
Set the nonlinear equation system that will be solved.
PETScNonlinearSolver(GlobalLinearSolver &linear_solver, int maxiter, std::string prefix)
bool _compensate_non_equilibrium_initial_residuum
std::size_t _jacobian_id
ID of the Jacobian matrix.
void compensateNonEquilibriumInitialResiduum(bool const value)
std::size_t _residual_id
ID of the residual vector.
void calculateNonEquilibriumInitialResiduum(std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id) override
virtual GlobalVector & getVector(std::size_t &id)=0
Get an uninitialized vector with the given id.
virtual void releaseVector(GlobalVector const &x)=0
@ FirstOrderImplicitQuasilinear
void copy(PETScVector const &x, PETScVector &y)
void finalizeVectorAssembly(VEC_T &)
General function to finalize the vector assembly.
PetscErrorCode updateResidual(SNES, Vec x, Vec petsc_r, void *petsc_context)
PetscErrorCode updateJacobian(SNES, Vec, Mat J, Mat, void *petsc_context)
static NUMLIB_EXPORT MatrixProvider & provider
static NUMLIB_EXPORT VectorProvider & provider
Status of the non-linear solver.
std::vector< GlobalVector * > & x
std::vector< GlobalVector * > const & x_prev