Skip to content

Commit

Permalink
refactor implicit update
Browse files Browse the repository at this point in the history
  • Loading branch information
@landinjm
landinjm committed Oct 18, 2024
1 parent 3d7e410 commit d235d9b
Show file tree
Hide file tree
Showing 3 changed files with 260 additions and 523 deletions.
274 changes: 4 additions & 270 deletions applications/allenCahn_conserved/customPDE.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -211,277 +211,11 @@ customPDE<dim, degree>::solveIncrement(bool skip_time_dependent)
this->pcout << buffer;
}

LinearAlgebra::distributed::Vector<double> solution_diff =
*this->solutionSet[fieldIndex];

// apply Dirichlet BC's
// This clears the residual where we want to apply Dirichlet
// BCs, otherwise the solver sees a positive residual
this->constraintsDirichletSet[fieldIndex]->set_zero(
*this->residualSet[fieldIndex]);

// solver controls
double tol_value;
if (MatrixFreePDE<dim, degree>::userInputs.linear_solver_parameters
.getToleranceType(fieldIndex) == ABSOLUTE_RESIDUAL)
{
tol_value =
MatrixFreePDE<dim, degree>::userInputs.linear_solver_parameters
.getToleranceValue(fieldIndex);
}
else
{
tol_value =
MatrixFreePDE<dim, degree>::userInputs.linear_solver_parameters
.getToleranceValue(fieldIndex) *
this->residualSet[fieldIndex]->l2_norm();
}

SolverControl solver_control(
MatrixFreePDE<dim, degree>::userInputs.linear_solver_parameters
.getMaxIterations(fieldIndex),
tol_value);

// Currently the only allowed solver is SolverCG, the
// SolverType input variable is a dummy
SolverCG<vectorType> solver(solver_control);

// solve
try
{
if (this->fields[fieldIndex].type == SCALAR)
{
this->dU_scalar = 0.0;
solver.solve(*this,
this->dU_scalar,
*this->residualSet[fieldIndex],
IdentityMatrix(
this->solutionSet[fieldIndex]->size()));
}
else
{
this->dU_vector = 0.0;
solver.solve(*this,
this->dU_vector,
*this->residualSet[fieldIndex],
IdentityMatrix(
this->solutionSet[fieldIndex]->size()));
}
}
catch (...)
{
this->pcout << "\nWarning: linear solver did not converge as per "
"set tolerances. consider increasing the maximum "
"number of iterations or decreasing the solver "
"tolerance.\n";
}

if (userInputs.var_nonlinear[fieldIndex])
{
// Now that we have the calculated change in the solution,
// we need to select a damping coefficient
double damping_coefficient;

if (MatrixFreePDE<dim, degree>::userInputs
.nonlinear_solver_parameters.getBacktrackDampingFlag(
fieldIndex))
{
vectorType solutionSet_old = *this->solutionSet[fieldIndex];
double residual_old = this->residualSet[fieldIndex]->l2_norm();

damping_coefficient = 1.0;
bool damping_coefficient_found = false;
while (!damping_coefficient_found)
{
if (this->fields[fieldIndex].type == SCALAR)
{
this->solutionSet[fieldIndex]->sadd(1.0,
damping_coefficient,
this->dU_scalar);
}
else
{
this->solutionSet[fieldIndex]->sadd(1.0,
damping_coefficient,
this->dU_vector);
}

this->computeNonexplicitRHS();

for (const auto &it : *this->valuesDirichletSet[fieldIndex])
{
if (this->residualSet[fieldIndex]->in_local_range(
it.first))
{
(*this->residualSet[fieldIndex])(it.first) = 0.0;
}
}

double residual_new =
this->residualSet[fieldIndex]->l2_norm();

if (this->currentIncrement % userInputs.skip_print_steps ==
0)
{
this->pcout << " Old residual: " << residual_old
<< " Damping Coeff: " << damping_coefficient
<< " New Residual: " << residual_new
<< std::endl;
}

// An improved approach would use the
// Armijo–Goldstein condition to ensure a
// sufficent decrease in the residual. This way is
// just scales the residual.
if ((residual_new <
(residual_old *
MatrixFreePDE<dim, degree>::userInputs
.nonlinear_solver_parameters
.getBacktrackResidualDecreaseCoeff(fieldIndex))) ||
damping_coefficient < 1.0e-4)
{
damping_coefficient_found = true;
}
else
{
damping_coefficient *=
MatrixFreePDE<dim, degree>::userInputs
.nonlinear_solver_parameters
.getBacktrackStepModifier(fieldIndex);
*this->solutionSet[fieldIndex] = solutionSet_old;
}
}
}
else
{
damping_coefficient =
MatrixFreePDE<dim, degree>::userInputs
.nonlinear_solver_parameters.getDefaultDampingCoefficient(
fieldIndex);

if (this->fields[fieldIndex].type == SCALAR)
{
this->solutionSet[fieldIndex]->sadd(1.0,
damping_coefficient,
this->dU_scalar);
}
else
{
this->solutionSet[fieldIndex]->sadd(1.0,
damping_coefficient,
this->dU_vector);
}
}

if (this->currentIncrement % userInputs.skip_print_steps == 0)
{
double dU_norm;
if (this->fields[fieldIndex].type == SCALAR)
{
dU_norm = this->dU_scalar.l2_norm();
}
else
{
dU_norm = this->dU_vector.l2_norm();
}
snprintf(buffer,
sizeof(buffer),
"field '%2s' [linear solve]: initial "
"residual:%12.6e, current residual:%12.6e, "
"nsteps:%u, tolerance criterion:%12.6e, "
"solution: %12.6e, dU: %12.6e\n",
this->fields[fieldIndex].name.c_str(),
this->residualSet[fieldIndex]->l2_norm(),
solver_control.last_value(),
solver_control.last_step(),
solver_control.tolerance(),
this->solutionSet[fieldIndex]->l2_norm(),
dU_norm);
this->pcout << buffer;
}

// Check to see if this individual variable has converged
if (MatrixFreePDE<dim, degree>::userInputs
.nonlinear_solver_parameters.getToleranceType(fieldIndex) ==
ABSOLUTE_SOLUTION_CHANGE)
{
double diff;

if (this->fields[fieldIndex].type == SCALAR)
{
diff = this->dU_scalar.l2_norm();
}
else
{
diff = this->dU_vector.l2_norm();
}
if (this->currentIncrement % userInputs.skip_print_steps == 0)
{
this->pcout << "Relative difference between "
"nonlinear iterations: "
<< diff << " " << nonlinear_it_index << " "
<< this->currentIncrement << std::endl;
}
nonlinear_it_converged =
this->updateImplicitSolution(fieldIndex, nonlinear_it_index);

if (diff > MatrixFreePDE<dim, degree>::userInputs
.nonlinear_solver_parameters.getToleranceValue(
fieldIndex) &&
nonlinear_it_index <
MatrixFreePDE<dim, degree>::userInputs
.nonlinear_solver_parameters.getMaxIterations())
{
nonlinear_it_converged = false;
}
}
else
{
std::cerr << "PRISMS-PF Error: Nonlinear solver tolerance "
"types other than ABSOLUTE_CHANGE have yet to "
"be implemented."
<< std::endl;
}
}
else
{
if (nonlinear_it_index == 0)
{
if (this->fields[fieldIndex].type == SCALAR)
{
*this->solutionSet[fieldIndex] += this->dU_scalar;
}
else
{
*this->solutionSet[fieldIndex] += this->dU_vector;
}

if (this->currentIncrement % userInputs.skip_print_steps == 0)
{
double dU_norm;
if (this->fields[fieldIndex].type == SCALAR)
{
dU_norm = this->dU_scalar.l2_norm();
}
else
{
dU_norm = this->dU_vector.l2_norm();
}
snprintf(buffer,
sizeof(buffer),
"field '%2s' [linear solve]: initial "
"residual:%12.6e, current residual:%12.6e, "
"nsteps:%u, tolerance criterion:%12.6e, "
"solution: %12.6e, dU: %12.6e\n",
this->fields[fieldIndex].name.c_str(),
this->residualSet[fieldIndex]->l2_norm(),
solver_control.last_value(),
solver_control.last_step(),
solver_control.tolerance(),
this->solutionSet[fieldIndex]->l2_norm(),
dU_norm);
this->pcout << buffer;
}
}
}
// Apply Boundary conditions
this->applyBCs(fieldIndex);
}
else if (this->fields[fieldIndex].pdetype == AUXILIARY)
{
Expand Down
4 changes: 4 additions & 0 deletions include/matrixFreePDE.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -262,6 +262,10 @@ class MatrixFreePDE : public Subscriptor
void
updateExplicitSolution(unsigned int fieldIndex);

/*Method to compute an implicit timestep*/
bool
updateImplicitSolution(unsigned int fieldIndex, unsigned int nonlinear_it_index);

/*Method to apply boundary conditions*/
void
applyBCs(unsigned int fieldIndex);
Expand Down
Loading

0 comments on commit d235d9b

Please sign in to comment.