From e427dd7e66483409eec8f8f06e909be916eb7c5b Mon Sep 17 00:00:00 2001 From: Xander <102053371+fractalsbyx@users.noreply.github.com> Date: Thu, 16 Jan 2025 12:49:38 -0500 Subject: [PATCH] reorganize solveIncrement Broke solveIncrement into four smaller sections: 1. Aux vars that are not nonlinear (sequential) 2. Explicit vars (concurrent) 3. Linear solves (sequential) 4. Nonlinear solve (concurrent) TODO/NOTES: I broke the skip_time_dependent thing. Changes not yet propogated. Need some documentation. Prefer to pass var_index into solver rather than use member variable currentFieldIndex. Maybe removed need for member variable currentIncrement. computeExplicitRHS-like functions are commented with "HERE" to be replaced --- include/core/matrixFreePDE.h | 64 ++- src/core/solvers/solveIncrement.cc | 855 ++++++++++++++--------------- 2 files changed, 464 insertions(+), 455 deletions(-) diff --git a/include/core/matrixFreePDE.h b/include/core/matrixFreePDE.h index 07724f42..75c70940 100644 --- a/include/core/matrixFreePDE.h +++ b/include/core/matrixFreePDE.h @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -32,6 +33,7 @@ #include #include #include +#include #include #include @@ -171,6 +173,10 @@ class MatrixFreePDE : public Subscriptor std::vector> simplified_grain_representations; + // =========================================================================== + // METHODS FOR SOLVING + // =========================================================================== + /** * Method to solve each time increment of a time-dependent problem. For * time-independent problems this method is called only once. This method @@ -179,7 +185,51 @@ class MatrixFreePDE : public Subscriptor * problems, Implicit (matrix-free) solver for Elliptic problems. */ virtual void - solveIncrement(bool skip_time_dependent); + solveIncrement(uint current_increment, bool skip_time_dependent); + + void + linearSolveOnce(unsigned int var_index, SolverControl &solver_control); + void + auxiliaryOnce(unsigned int var_index, uint current_increment); + bool + nonlinearIncrement(unsigned int var_index, + SolverControl &solver_control, + uint current_increment); + bool + auxiliaryIncrement(unsigned int var_index, uint current_increment); + void + print_explicit_update(uint var_index, uint current_increment); + void + print_linear_update(uint var_index, + uint current_increment, + SolverControl &solver_control); + void + print_nonlinear_update(uint var_index, uint current_increment); + void + print_nonlinear_status(uint var_index, + uint current_increment, + uint nonlinear_iteration_index, + double diff); + + /*Method to compute an explicit timestep*/ + void + updateExplicitSolution(unsigned int var_index); + + void + updateLinearSolution(unsigned int var_index); + + /*Method to compute an implicit timestep*/ + bool + updateImplicitSolution(unsigned int fieldIndex, unsigned int nonlinear_iteration_index); + + /*Method to apply boundary conditions*/ + void + applyBCs(unsigned int var_index); + + // =========================================================================== + // + // =========================================================================== + /* Method to write solution fields to vtu and pvtu (parallel) files. * * This method can be enabled/disabled by setting the flag writeOutput to @@ -263,18 +313,6 @@ class MatrixFreePDE : public Subscriptor void computeInvM(); - /*Method to compute an explicit timestep*/ - void - updateExplicitSolution(unsigned int fieldIndex); - - /*Method to compute an implicit timestep*/ - bool - updateImplicitSolution(unsigned int fieldIndex, unsigned int nonlinear_iteration_index); - - /*Method to apply boundary conditions*/ - void - applyBCs(unsigned int fieldIndex); - /** * \brief Compute element volume for the triangulation */ diff --git a/src/core/solvers/solveIncrement.cc b/src/core/solvers/solveIncrement.cc index afdd1f64..a6564324 100644 --- a/src/core/solvers/solveIncrement.cc +++ b/src/core/solvers/solveIncrement.cc @@ -1,218 +1,354 @@ #include +#include + +#include "core/typeEnums.h" #include #include #include #include +#include // solve each time increment template void -MatrixFreePDE::solveIncrement(bool skip_time_dependent) +MatrixFreePDE::solveIncrement(uint current_increment, + bool skip_time_dependent) { + // log time + computing_timer.enter_subsection("matrixFreePDE: solveIncrements"); char buffer[200]; - // Get the RHS of the explicit equations - if (hasExplicitEquation && !skip_time_dependent) + // Solve and update auxiliary fields that are not nonlinear (sequentially) + for (const auto &[var_index, variable] : var_attributes) { - // computeExplicitRHS(); + if (variable.eq_type == AUXILIARY && !variable.is_nonlinear) + { + auxiliaryOnce(var_index); + } } - // solve for each field - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + // Solve and update explicit fields (concurrently) + // if (hasExplicitEquation && !skip_time_dependent) + computeExplicitRHS(); // HERE + for (const auto &[var_index, variable] : var_attributes) { - currentFieldIndex = fieldIndex; // Used in computeLHS() + if (variable.eq_type == EXPLICIT_TIME_DEPENDENT) + { + updateExplicitSolution(var_index); + // Apply Boundary conditions + applyBCs(var_index); + print_explicit_update(var_index, current_increment); + } + } - // Parabolic (first order derivatives in time) fields - if (fields[fieldIndex].pdetype == EXPLICIT_TIME_DEPENDENT && !skip_time_dependent) + // Solve and update linear-equation fields (sequentially) + for (const auto &[var_index, variable] : var_attributes) + { + if (variable.eq_type != EXPLICIT_TIME_DEPENDENT && !variable.is_nonlinear) { - updateExplicitSolution(fieldIndex); + SolverControl solver_control; + compute_nonexplicit_RHS(var_index); // HERE + linearSolveOnce(var_index, solver_control); // lhs + updateLinearSolution(var_index, solver_control); + print_linear_update(var_index, solver_control); // Apply Boundary conditions - applyBCs(fieldIndex); + applyBCs(var_index); + print_update(var_index, current_increment); + } + } - // Print update to screen and confirm that solution isn't nan - if (currentIncrement % userInputs.skip_print_steps == 0) + // Solve and update nonlinear-equation fields (concurrently) + bool nonlinear_converged = false; + uint nonlinear_iteration_index = 0; + uint max_nonlinear_iterations = + userInputs.nonlinear_solver_parameters.getMaxIterations(); + SolverControl solver_control; + while (!nonlinear_converged && nonlinear_iteration_index < max_nonlinear_iterations) + { + nonlinear_converged = true; + for (const auto &[var_index, variable] : var_attributes) + { + if (!variable.is_nonlinear) + { + continue; + } + if (variable.eq_type == TIME_INDEPENDENT) + { + print_nonlinear_update(var_index, current_increment); + if (!nonlinearIncrement(var_index, solver_control)) + { + nonlinear_converged = false; + } + // Apply Boundary conditions (could be placed in nonlinearIncrement) + applyBCs(var_index); + } + else if (variable.eq_type == AUXILIARY) + { + if (!auxiliaryIncrement(var_index, solver_control)) + { + nonlinear_converged = false; + } + // Apply Boundary conditions (could be placed in auxiliaryIncrement) + applyBCs(var_index); + } + // check if solution is nan + if (!numbers::is_finite(solutionSet[var_index]->l2_norm())) { - double solution_L2_norm = solutionSet[fieldIndex]->l2_norm(); - snprintf(buffer, sizeof(buffer), - "field '%2s' [explicit solve]: current solution: " - "%12.6e, current residual:%12.6e\n", - fields[fieldIndex].name.c_str(), - solution_L2_norm, - residualSet[fieldIndex]->l2_norm()); + "ERROR: field '%s' solution is NAN. exiting.\n\n", + fields[var_index].name.c_str()); conditionalOStreams::pout_base << buffer; - - if (!numbers::is_finite(solution_L2_norm)) - { - snprintf(buffer, - sizeof(buffer), - "ERROR: field '%s' solution is NAN. exiting.\n\n", - fields[fieldIndex].name.c_str()); - conditionalOStreams::pout_base << buffer; - exit(-1); - } + exit(-1); } } + nonlinear_iteration_index++; + } + if (!nonlinear_converged) + { + conditionalOStreams::pout_base << "\nWarning: nonlinear solver did not converge as " + "per set tolerances. consider increasing the " + "maximum number of iterations or decreasing the " + "solver tolerance.\n"; } + // log time + computing_timer.leave_subsection("matrixFreePDE: solveIncrements"); +} - // Now, update the non-explicit variables - // For the time being, this is just the elliptic equations, but implicit - // parabolic and auxilary equations should also be here - if (hasNonExplicitEquation) +template +void +MatrixFreePDE::linearSolveOnce(unsigned int var_index, + SolverControl &solver_control) +{ + // Apply Dirichlet BC's. This clears the residual where we want to apply Dirichlet BCs, + // otherwise the solver sees a positive residual + constraintsDirichletSet[var_index]->set_zero(*residualSet[var_index]); + + // Grab solver controls + double tol_value = NAN; + if (userInputs.linear_solver_parameters.getToleranceType(var_index) == + ABSOLUTE_RESIDUAL) { - bool nonlinear_iteration_converged = false; - unsigned int nonlinear_iteration_index = 0; + tol_value = userInputs.linear_solver_parameters.getToleranceValue(var_index); + } + else + { + tol_value = userInputs.linear_solver_parameters.getToleranceValue(var_index) * + residualSet[var_index]->l2_norm(); + } + solver_control = + SolverControl(userInputs.linear_solver_parameters.getMaxIterations(var_index), + tol_value); + SolverCG> solver(solver_control); - while (!nonlinear_iteration_converged) + // Solve the linear system + try + { + if (fields[var_index].type == SCALAR) + { + dU_scalar = 0.0; + solver.solve(*this, + dU_scalar, + *residualSet[var_index], + IdentityMatrix(solutionSet[var_index]->size())); + } + else { - nonlinear_iteration_converged = true; + dU_vector = 0.0; + solver.solve(*this, + dU_vector, + *residualSet[var_index], + IdentityMatrix(solutionSet[var_index]->size())); + } + } + catch (...) + { + conditionalOStreams::pout_base << "\nWarning: linear solver did not converge as " + "per set tolerances. consider increasing the " + "maximum number of iterations or decreasing the " + "solver tolerance.\n"; + } +} + +template +void +MatrixFreePDE::auxiliaryOnce(unsigned int var_index, uint current_increment) +{ + computeNonexplicitRHS(var_index); // HERE + updateExplicitSolution(var_index); - // Update residualSet for the non-explicitly updated variables - // computeNonexplicitRHS(); + // Apply Boundary conditions + applyBCs(var_index); - for (const auto &[fieldIndex, variable] : var_attributes) + // Print update to screen + if (current_increment % userInputs.skip_print_steps == 0) + { + char buffer[200]; + snprintf(buffer, + sizeof(buffer), + "field '%2s' [auxiliary solve]: current solution: " + "%12.6e, current residual:%12.6e\n", + fields[var_index].name.c_str(), + solutionSet[var_index]->l2_norm(), + residualSet[var_index]->l2_norm()); + conditionalOStreams::pout_base << buffer; + } +} + +// Nonlinear increment for time-independent fields +template +bool +MatrixFreePDE::nonlinearIncrement(unsigned int var_index, + SolverControl &solver_control, + uint current_increment) +{ + linearSolveOnce(var_index, solver_control); // lhs + + uint nonlinear_iteration_index = 0; + double damping_coefficient = NAN; + if (userInputs.nonlinear_solver_parameters.getBacktrackDampingFlag(var_index)) + { + dealii::LinearAlgebra::distributed::Vector solutionSet_old = + *solutionSet[var_index]; + double residual_old = residualSet[var_index]->l2_norm(); + + damping_coefficient = 1.0; + bool damping_coefficient_found = false; + while (!damping_coefficient_found) + { + if (fields[var_index].type == SCALAR) { - currentFieldIndex = fieldIndex; // Used in computeLHS() + solutionSet[var_index]->sadd(1.0, damping_coefficient, dU_scalar); + } + else + { + solutionSet[var_index]->sadd(1.0, damping_coefficient, dU_vector); + } - if ((fields[fieldIndex].pdetype == IMPLICIT_TIME_DEPENDENT && - !skip_time_dependent) || - fields[fieldIndex].pdetype == TIME_INDEPENDENT) - { - if (currentIncrement % userInputs.skip_print_steps == 0 && - variable.is_nonlinear) - { - snprintf(buffer, - sizeof(buffer), - "field '%2s' [nonlinear solve]: current " - "solution: %12.6e, current residual:%12.6e\n", - fields[fieldIndex].name.c_str(), - solutionSet[fieldIndex]->l2_norm(), - residualSet[fieldIndex]->l2_norm()); - conditionalOStreams::pout_base << buffer; - } - - nonlinear_iteration_converged = - updateImplicitSolution(fieldIndex, nonlinear_iteration_index); - - // Apply Boundary conditions - applyBCs(fieldIndex); - } - else if (fields[fieldIndex].pdetype == AUXILIARY) - { - if (variable.is_nonlinear || nonlinear_iteration_index == 0) - { - // If the equation for this field is nonlinear, save the old - // solution - if (variable.is_nonlinear) - { - if (fields[fieldIndex].type == SCALAR) - { - dU_scalar = *solutionSet[fieldIndex]; - } - else - { - dU_vector = *solutionSet[fieldIndex]; - } - } - - updateExplicitSolution(fieldIndex); - - // Apply Boundary conditions - applyBCs(fieldIndex); - - // Print update to screen - if (currentIncrement % userInputs.skip_print_steps == 0) - { - snprintf(buffer, - sizeof(buffer), - "field '%2s' [auxiliary solve]: current solution: " - "%12.6e, current residual:%12.6e\n", - fields[fieldIndex].name.c_str(), - solutionSet[fieldIndex]->l2_norm(), - residualSet[fieldIndex]->l2_norm()); - conditionalOStreams::pout_base << buffer; - } - - // Check to see if this individual variable has converged - if (variable.is_nonlinear) - { - if (userInputs.nonlinear_solver_parameters.getToleranceType( - fieldIndex) == ABSOLUTE_SOLUTION_CHANGE) - { - double diff = NAN; - - if (fields[fieldIndex].type == SCALAR) - { - dU_scalar -= *solutionSet[fieldIndex]; - diff = dU_scalar.l2_norm(); - } - else - { - dU_vector -= *solutionSet[fieldIndex]; - diff = dU_vector.l2_norm(); - } - if (currentIncrement % userInputs.skip_print_steps == 0) - { - snprintf(buffer, - sizeof(buffer), - " field '%2s' [nonlinear solve] current " - "increment: %u, nonlinear " - "iteration: " - "%u, dU: %12.6e\n", - fields[fieldIndex].name.c_str(), - currentIncrement, - nonlinear_iteration_index, - diff); - conditionalOStreams::pout_base << buffer; - } - - if (diff > userInputs.nonlinear_solver_parameters - .getToleranceValue(fieldIndex) && - nonlinear_iteration_index < - userInputs.nonlinear_solver_parameters - .getMaxIterations()) - { - nonlinear_iteration_converged = false; - } - } - else - { - AssertThrow( - false, - FeatureNotImplemented( - "Nonlinear solver tolerances besides ABSOLUTE_CHANGE")); - } - } - } - } + computeNonexplicitRHS(); // HERE - // check if solution is nan - if (!numbers::is_finite(solutionSet[fieldIndex]->l2_norm())) + for (const auto &it : *valuesDirichletSet[var_index]) + { + if (residualSet[var_index]->in_local_range(it.first)) { - snprintf(buffer, - sizeof(buffer), - "ERROR: field '%s' solution is NAN. exiting.\n\n", - fields[fieldIndex].name.c_str()); - conditionalOStreams::pout_base << buffer; - exit(-1); + (*residualSet[var_index])(it.first) = 0.0; } } - nonlinear_iteration_index++; + double residual_new = residualSet[var_index]->l2_norm(); + + if (current_increment % userInputs.skip_print_steps == 0) + { + conditionalOStreams::pout_base << "\tOld residual: " << residual_old + << " Damping Coeff: " << damping_coefficient + << " New Residual: " << residual_new << "\n"; + } + + // 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 * + userInputs.nonlinear_solver_parameters.getBacktrackResidualDecreaseCoeff( + var_index))) || + damping_coefficient < 1.0e-4) + { + damping_coefficient_found = true; + } + else + { + damping_coefficient *= + userInputs.nonlinear_solver_parameters.getBacktrackStepModifier( + var_index); + *solutionSet[var_index] = solutionSet_old; + } } } + else + { + damping_coefficient = + userInputs.nonlinear_solver_parameters.getDefaultDampingCoefficient(var_index); + + if (fields[var_index].type == SCALAR) + { + solutionSet[var_index]->sadd(1.0, damping_coefficient, dU_scalar); + } + else + { + solutionSet[var_index]->sadd(1.0, damping_coefficient, dU_vector); + } + } + + // Print linear???? + print_linear_update(var_index, solver_control); + + // Check to see if this individual variable has converged + bool nonlinear_converged = false; + AssertThrow(userInputs.nonlinear_solver_parameters.getToleranceType(var_index) == + ABSOLUTE_SOLUTION_CHANGE, + FeatureNotImplemented( + "Nonlinear solver tolerances besides ABSOLUTE_CHANGE")); + double diff = NAN; + if (fields[var_index].type == SCALAR) + { + diff = dU_scalar.l2_norm(); + } + else + { + diff = dU_vector.l2_norm(); + } + print_nonlinear_status(var_index, current_increment, nonlinear_iteration_index, diff); + return !(diff > userInputs.nonlinear_solver_parameters.getToleranceValue(var_index)); +} + +// Nonlinear increment for auxiliary fields +template +bool +MatrixFreePDE::auxiliaryIncrement(unsigned int var_index, + uint current_increment) +{ + // Save the old solution + if (fields[var_index].type == SCALAR) + { + dU_scalar = *solutionSet[var_index]; + } + else + { + dU_vector = *solutionSet[var_index]; + } + + auxiliaryOnce(var_index); + + // Check to see if this individual variable has converged + AssertThrow(userInputs.nonlinear_solver_parameters.getToleranceType(var_index) == + ABSOLUTE_SOLUTION_CHANGE, + FeatureNotImplemented( + "Nonlinear solver tolerances besides ABSOLUTE_CHANGE")); + double diff = NAN; + if (fields[var_index].type == SCALAR) + { + dU_scalar -= *solutionSet[var_index]; + diff = dU_scalar.l2_norm(); + } + else + { + dU_vector -= *solutionSet[var_index]; + diff = dU_vector.l2_norm(); + } + print_nonlinear_update(var_index, current_increment); + return !(diff > userInputs.nonlinear_solver_parameters.getToleranceValue(var_index)); } // Application of boundary conditions template void -MatrixFreePDE::applyBCs(unsigned int fieldIndex) +MatrixFreePDE::applyBCs(unsigned int var_index) { // Add Neumann BCs - if (fields[fieldIndex].hasNeumannBCs) + if (fields[var_index].hasNeumannBCs) { // Currently commented out because it isn't working yet // applyNeumannBCs(); @@ -220,344 +356,179 @@ MatrixFreePDE::applyBCs(unsigned int fieldIndex) // Set the Dirichelet values (hanging node constraints don't need to be distributed // every time step, only at output) - if (fields[fieldIndex].hasDirichletBCs) + if (fields[var_index].hasDirichletBCs) { // Apply non-uniform Dirlichlet_BCs to the current field - if (fields[fieldIndex].hasnonuniformDirichletBCs) + if (fields[var_index].hasnonuniformDirichletBCs) { - DoFHandler *dof_handler = nullptr; - dof_handler = dofHandlersSet_nonconst.at(currentFieldIndex); - IndexSet *locally_relevant_dofs = nullptr; - locally_relevant_dofs = locally_relevant_dofsSet_nonconst.at(currentFieldIndex); + DoFHandler *dof_handler = dofHandlersSet_nonconst.at(currentvar_index); + IndexSet *locally_relevant_dofs = + locally_relevant_dofsSet_nonconst.at(currentvar_index); locally_relevant_dofs->clear(); DoFTools::extract_locally_relevant_dofs(*dof_handler, *locally_relevant_dofs); - AffineConstraints *constraintsDirichlet = nullptr; - constraintsDirichlet = constraintsDirichletSet_nonconst.at(currentFieldIndex); + AffineConstraints *constraintsDirichlet = + constraintsDirichletSet_nonconst.at(currentvar_index); constraintsDirichlet->clear(); constraintsDirichlet->reinit(*locally_relevant_dofs); applyDirichletBCs(); constraintsDirichlet->close(); } // Distribute for Uniform or Non-Uniform Dirichlet BCs - constraintsDirichletSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); + constraintsDirichletSet[var_index]->distribute(*solutionSet[var_index]); } - solutionSet[fieldIndex]->update_ghost_values(); + solutionSet[var_index]->update_ghost_values(); } // Explicit time step for matrixfree solve template void -MatrixFreePDE::updateExplicitSolution(unsigned int fieldIndex) +MatrixFreePDE::updateExplicitSolution(unsigned int var_index) { // Explicit-time step each DOF // Takes advantage of knowledge that the length of solutionSet and residualSet // is an integer multiple of the length of invM for vector variables - if (fields[fieldIndex].type == SCALAR) + if (fields[var_index].type == SCALAR) { unsigned int invM_size = invMscalar.locally_owned_size(); - for (unsigned int dof = 0; dof < solutionSet[fieldIndex]->locally_owned_size(); + for (unsigned int dof = 0; dof < solutionSet[var_index]->locally_owned_size(); ++dof) { - solutionSet[fieldIndex]->local_element(dof) = + solutionSet[var_index]->local_element(dof) = invMscalar.local_element(dof % invM_size) * - residualSet[fieldIndex]->local_element(dof); + residualSet[var_index]->local_element(dof); } } - else if (fields[fieldIndex].type == VECTOR) + else if (fields[var_index].type == VECTOR) { unsigned int invM_size = invMvector.locally_owned_size(); - for (unsigned int dof = 0; dof < solutionSet[fieldIndex]->locally_owned_size(); + for (unsigned int dof = 0; dof < solutionSet[var_index]->locally_owned_size(); ++dof) { - solutionSet[fieldIndex]->local_element(dof) = + solutionSet[var_index]->local_element(dof) = invMvector.local_element(dof % invM_size) * - residualSet[fieldIndex]->local_element(dof); + residualSet[var_index]->local_element(dof); } } } template -bool -MatrixFreePDE::updateImplicitSolution(unsigned int fieldIndex, - unsigned int nonlinear_iteration_index) +void +MatrixFreePDE::updateLinearSolution(unsigned int var_index) { - char buffer[200]; - - // Assume convergence criterion is met, unless otherwise proven later on. - bool nonlinear_iteration_converged = true; - - // Apply Dirichlet BC's. This clears the residual where we want to apply Dirichlet BCs, - // otherwise the solver sees a positive residual - constraintsDirichletSet[fieldIndex]->set_zero(*residualSet[fieldIndex]); - - // Grab solver controls - double tol_value = NAN; - if (userInputs.linear_solver_parameters.getToleranceType(fieldIndex) == - ABSOLUTE_RESIDUAL) + if (fields[var_index].type == SCALAR) { - tol_value = userInputs.linear_solver_parameters.getToleranceValue(fieldIndex); + *solutionSet[var_index] += dU_scalar; } else { - tol_value = userInputs.linear_solver_parameters.getToleranceValue(fieldIndex) * - residualSet[fieldIndex]->l2_norm(); + *solutionSet[var_index] += dU_vector; } +} - SolverControl solver_control(userInputs.linear_solver_parameters.getMaxIterations( - fieldIndex), - tol_value); - - // Currently the only allowed solver is SolverCG, the - // SolverType input variable is a dummy - SolverCG> solver(solver_control); - - // Solve - try - { - if (fields[fieldIndex].type == SCALAR) - { - dU_scalar = 0.0; - // solver.solve(*this, - // dU_scalar, - // *residualSet[fieldIndex], - // IdentityMatrix(solutionSet[fieldIndex]->size())); - } - else - { - dU_vector = 0.0; - // solver.solve(*this, - // dU_vector, - // *residualSet[fieldIndex], - // IdentityMatrix(solutionSet[fieldIndex]->size())); - } - } - catch (...) - { - conditionalOStreams::pout_base << "\nWarning: linear solver did not converge as " - "per set tolerances. consider increasing the " - "maximum number of iterations or decreasing the " - "solver tolerance.\n"; - } +template +void +MatrixFreePDE::print_explicit_update(uint var_index, uint current_increment) +{ + char buffer[200]; - if (var_attributes.at(fieldIndex).is_nonlinear) + if (current_increment % userInputs.skip_print_steps == 0) { - // Now that we have the calculated change in the solution, - // we need to select a damping coefficient - double damping_coefficient = NAN; - - if (userInputs.nonlinear_solver_parameters.getBacktrackDampingFlag(fieldIndex)) - { - dealii::LinearAlgebra::distributed::Vector solutionSet_old = - *solutionSet[fieldIndex]; - double residual_old = residualSet[fieldIndex]->l2_norm(); - - damping_coefficient = 1.0; - bool damping_coefficient_found = false; - while (!damping_coefficient_found) - { - if (fields[fieldIndex].type == SCALAR) - { - solutionSet[fieldIndex]->sadd(1.0, damping_coefficient, dU_scalar); - } - else - { - solutionSet[fieldIndex]->sadd(1.0, damping_coefficient, dU_vector); - } - - // computeNonexplicitRHS(); - - for (const auto &it : *valuesDirichletSet[fieldIndex]) - { - if (residualSet[fieldIndex]->in_local_range(it.first)) - { - (*residualSet[fieldIndex])(it.first) = 0.0; - } - } - - double residual_new = residualSet[fieldIndex]->l2_norm(); - - if (currentIncrement % userInputs.skip_print_steps == 0) - { - conditionalOStreams::pout_base - << " Old residual: " << residual_old - << " Damping Coeff: " << damping_coefficient - << " New Residual: " << residual_new << "\n"; - } - - // 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 * userInputs.nonlinear_solver_parameters - .getBacktrackResidualDecreaseCoeff(fieldIndex))) || - damping_coefficient < 1.0e-4) - { - damping_coefficient_found = true; - } - else - { - damping_coefficient *= - userInputs.nonlinear_solver_parameters.getBacktrackStepModifier( - fieldIndex); - *solutionSet[fieldIndex] = solutionSet_old; - } - } - } - else - { - damping_coefficient = - userInputs.nonlinear_solver_parameters.getDefaultDampingCoefficient( - fieldIndex); - - if (fields[fieldIndex].type == SCALAR) - { - solutionSet[fieldIndex]->sadd(1.0, damping_coefficient, dU_scalar); - } - else - { - solutionSet[fieldIndex]->sadd(1.0, damping_coefficient, dU_vector); - } - } - - if (currentIncrement % userInputs.skip_print_steps == 0) + double solution_L2_norm = solutionSet[var_index]->l2_norm(); + + snprintf(buffer, + sizeof(buffer), + "field '%2s' [explicit solve]: current solution: " + "%12.6e, current residual:%12.6e\n", + fields[var_index].name.c_str(), + solution_L2_norm, + residualSet[var_index]->l2_norm()); + conditionalOStreams::pout_base << buffer; + + if (!numbers::is_finite(solution_L2_norm)) { - double dU_norm = NAN; - if (fields[fieldIndex].type == SCALAR) - { - dU_norm = dU_scalar.l2_norm(); - } - else - { - dU_norm = 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", - fields[fieldIndex].name.c_str(), - residualSet[fieldIndex]->l2_norm(), - solver_control.last_value(), - solver_control.last_step(), - solver_control.tolerance(), - solutionSet[fieldIndex]->l2_norm(), - dU_norm); + "ERROR: field '%s' solution is NAN. exiting.\n\n", + fields[var_index].name.c_str()); conditionalOStreams::pout_base << buffer; + exit(-1); } + } +} - // Check to see if this individual variable has converged - if (userInputs.nonlinear_solver_parameters.getToleranceType(fieldIndex) == - ABSOLUTE_SOLUTION_CHANGE) +template +void +MatrixFreePDE::print_linear_update(uint var_index, + uint current_increment, + SolverControl &solver_control) +{ + if (current_increment % userInputs.skip_print_steps == 0) + { + double dU_norm = NAN; + if (fields[var_index].type == SCALAR) { - double diff = NAN; - - if (fields[fieldIndex].type == SCALAR) - { - diff = dU_scalar.l2_norm(); - } - else - { - diff = dU_vector.l2_norm(); - } - if (currentIncrement % userInputs.skip_print_steps == 0) - { - snprintf(buffer, - sizeof(buffer), - " field '%2s' [nonlinear solve] current increment: %u, nonlinear " - "iteration: " - "%u, dU: %12.6e\n", - fields[fieldIndex].name.c_str(), - currentIncrement, - nonlinear_iteration_index, - diff); - conditionalOStreams::pout_base << buffer; - } - - if (diff > - userInputs.nonlinear_solver_parameters.getToleranceValue(fieldIndex) && - nonlinear_iteration_index < - userInputs.nonlinear_solver_parameters.getMaxIterations()) - { - nonlinear_iteration_converged = false; - } - else if (diff > - userInputs.nonlinear_solver_parameters.getToleranceValue(fieldIndex)) - { - conditionalOStreams::pout_base - << "\nWarning: nonlinear solver did not converge as " - "per set tolerances. consider increasing the " - "maximum number of iterations or decreasing the " - "solver tolerance.\n"; - } + dU_norm = dU_scalar.l2_norm(); } else { - AssertThrow(false, - FeatureNotImplemented( - "Nonlinear solver tolerances besides ABSOLUTE_CHANGE")); + dU_norm = dU_vector.l2_norm(); } + char buffer[200]; + 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", + fields[var_index].name.c_str(), + residualSet[var_index]->l2_norm(), + solver_control.last_value(), + solver_control.last_step(), + solver_control.tolerance(), + solutionSet[var_index]->l2_norm(), + dU_norm); + conditionalOStreams::pout_base << buffer; } - else - { - if (nonlinear_iteration_index == 0) - { - if (fields[fieldIndex].type == SCALAR) - { - *solutionSet[fieldIndex] += dU_scalar; - } - else - { - *solutionSet[fieldIndex] += dU_vector; - } - - if (currentIncrement % userInputs.skip_print_steps == 0) - { - double dU_norm = NAN; - if (fields[fieldIndex].type == SCALAR) - { - dU_norm = dU_scalar.l2_norm(); - } - else - { - dU_norm = 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", - fields[fieldIndex].name.c_str(), - residualSet[fieldIndex]->l2_norm(), - solver_control.last_value(), - solver_control.last_step(), - solver_control.tolerance(), - solutionSet[fieldIndex]->l2_norm(), - dU_norm); - conditionalOStreams::pout_base << buffer; - } - } - } - - return nonlinear_iteration_converged; } -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; +template +void +MatrixFreePDE::print_nonlinear_update(uint var_index, uint current_increment) +{ + char buffer[200]; -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; + if (!(current_increment % userInputs.skip_print_steps)) + { + snprintf(buffer, + sizeof(buffer), + "field '%2s' [nonlinear solve]: current " + "solution: %12.6e, current residual:%12.6e\n", + fields[var_index].name.c_str(), + solutionSet[var_index]->l2_norm(), + residualSet[var_index]->l2_norm()); + conditionalOStreams::pout_base << buffer; + } +} -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file +template +void +MatrixFreePDE::print_nonlinear_status(uint var_index, + uint current_increment, + uint nonlinear_iteration_index, + double diff) +{ + if (!(current_increment % userInputs.skip_print_steps)) + { + char buffer[200]; + snprintf(buffer, + sizeof(buffer), + " field '%2s' [nonlinear solve] current increment: %u, nonlinear " + "iteration: " + "%u, dU: %12.6e\n", + fields[var_index].name.c_str(), + current_increment, + nonlinear_iteration_index, + diff); + conditionalOStreams::pout_base << buffer; + } +} \ No newline at end of file