80 double rCurrentGuess[SIZE])
83 const double eps = 1e-6;
86 rCell.ComputeResidual(time, rCurrentGuess,
mResidual.data());
87 double norm_of_residual = norm_inf(
mResidual);
88 assert(!std::isnan(norm_of_residual));
89 double norm_of_update = 0.0;
93 rCell.ComputeJacobian(time, rCurrentGuess,
mJacobian);
99 norm_of_update = norm_inf(
mUpdate);
102 for (
unsigned i=0; i<SIZE; i++)
104 rCurrentGuess[i] -=
mUpdate[i];
106 double norm_of_previous_residual = norm_of_residual;
107 rCell.ComputeResidual(time, rCurrentGuess,
mResidual.data());
109 if (norm_of_residual > norm_of_previous_residual && norm_of_update > eps)
117 double relative_change_max = 0.0;
118 unsigned relative_change_direction = 0;
119 for (
unsigned i=0; i<SIZE; i++)
121 double relative_change = fabs(
mUpdate[i]/rCurrentGuess[i]);
122 if (relative_change > relative_change_max)
124 relative_change_max = relative_change;
125 relative_change_direction = i;
129 if (relative_change_max > 1.0)
132 rCurrentGuess[relative_change_direction] += 0.8*
mUpdate[relative_change_direction];
133 rCell.ComputeResidual(time, rCurrentGuess,
mResidual.data());
135 WARNING(
"Residual increasing and one direction changing radically - back tracking in that direction");
144 EXCEPTION(
"Newton method diverged in CardiacNewtonSolver::Solve()");
148 while (norm_of_update > eps);
152 if (norm_of_residual > 2e-10)
154 WARN_ONCE_ONLY(
"Newton iteration terminated because update vector norm is small, but residual norm is not small.");