Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor linear/nonlinear solves in solveIncrement.cc #265

Merged
merged 2 commits into from
Oct 21, 2024
Merged

Refactor linear/nonlinear solves in solveIncrement.cc #265

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d235d9b
refactor implicit update
landinjm Oct 18, 2024
0fc9366
changed bool default initialization to true/converged
landinjm Oct 21, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
274 changes: 4 additions & 270 deletions applications/allenCahn_conserved/customPDE.h
View file
Open in desktop
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
Open in desktop
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