Prevent first Newton iteration to be skipped (#236)

* Add final log for Newton
* Create output path for debug monitor
* Redirect call in debug monitor
* Revert time integration message. Unfortunately the message is before dt is changed.
* Add option to opt out early when hitting machine precision.

Author: termi-official

src/solver/logging.jl

@@ -19,7 +19,7 @@ function integration_step_monitor(integrator::SciMLBase.DEIntegrator, progress_m
     # (; id)       = progress_monitor
     (; tprev, t, dt, iter)    = integrator
     # push!(msgs, id => "$id: integrating on [$t, $(t+dt)] with Δt=$dt.")
-    @logmsg LogLevel(-100) "Integrating on [$tprev, $t]." iter=iter Δt=dt _group=:timeintegration
+    @logmsg LogLevel(-100) "Integrating on [$t, $(t+dt)]." iter=iter Δt=dt tprev=tprev _group=:timeintegration
 end
 
 function integration_finalize_monitor(integrator, progress_monitor::DefaultProgressMonitor)
@@ -39,15 +39,24 @@ function nonlinear_step_monitor(nlcache, time, f, u, progress_monitor::DefaultPr
     resnorm = norm(nlcache.residual)
     normΔu = norm(linear_solver_cache.u)
     if stats === nothing
-        @logmsg LogLevel(-100) "Nonlinear solve step" iter=iter resnorm=resnorm normΔu=normΔu _group=:nlsolve
+        @logmsg LogLevel(-100) "Nonlinear solve step" time=time iter=iter resnorm=resnorm normΔu=normΔu _group=:nlsolve
     else
-        @logmsg LogLevel(-100) "Nonlinear solve step" iter=iter stats=stats _group=:nlsolve
+        @logmsg LogLevel(-100) "Nonlinear solve step" time=time iter=iter resnorm=resnorm normΔu=normΔu stats=stats _group=:nlsolve
     end
 end
 
 function nonlinear_finalize_monitor(nlcache, time, f, progress_monitor::DefaultProgressMonitor)
     # (; id, msgs) = progress_monitor
-    # push!(msgs, id => "$id: done.")
+    (; iter, linear_solver_cache) = nlcache
+    stats = hasproperty(linear_solver_cache.cacheval, :stats) ? linear_solver_cache.cacheval.stats : nothing
+    # push!(msgs, id => "$id: $(nlcache.iter)\n\t||r||=$(norm(nlcache.residual)) ||Δu||=$(norm(linear_solver_cache.u))\n\t$stats")
+    resnorm = norm(nlcache.residual)
+    normΔu = norm(linear_solver_cache.u)
+    if stats === nothing
+        @logmsg LogLevel(-100) "Nonlinear solve converged" time=time iter=iter resnorm=resnorm normΔu=normΔu _group=:nlsolve
+    else
+        @logmsg LogLevel(-100) "Nonlinear solve converged" time=time iter=iter resnorm=resnorm normΔu=normΔu stats=stats _group=:nlsolve
+    end
 end
 
 #
@@ -71,12 +80,14 @@ end
 function nonlinear_step_monitor(cache, time, f, u, monitor::VTKNewtonMonitor)
     nonlinear_step_monitor(cache,time,f,u,monitor.inner_monitor)
 
+    mkpath(monitor.outdir)
     VTKGridFile(joinpath(monitor.outdir, "newton-monitor-t=$time-i=$(cache.iter).vtu"), f.dh) do vtk
         write_solution(vtk, f.dh, u)
         write_solution(vtk, f.dh, cache.linear_solver_cache.b, "_residual")
         write_solution(vtk, f.dh, cache.linear_solver_cache.u, "_increment")
     end
 end
 
-function nonlinear_finalize_monitor(nlcache, time, f, progress_monitor::VTKNewtonMonitor)
+function nonlinear_finalize_monitor(nlcache, time, f, monitor::VTKNewtonMonitor)
+    nonlinear_finalize_monitor(nlcache, time,f, monitor.inner_monitor)
 end

src/solver/nonlinear/newton_raphson.jl

@@ -68,7 +68,7 @@ function setup_solver_cache(f::AbstractSemidiscreteBlockedFunction, solver::Newt
     NewtonRaphsonSolverCache(op, residual, solver, inner_cache, T[], 0)
 end
 
-function nlsolve!(u::AbstractVector, f::AbstractSemidiscreteFunction, cache::NewtonRaphsonSolverCache, t)
+function nlsolve!(u::AbstractVector{T}, f::AbstractSemidiscreteFunction, cache::NewtonRaphsonSolverCache, t) where T
     @unpack op, residual, linear_solver_cache, Θks = cache
     monitor = cache.parameters.monitor
     cache.iter = -1
@@ -84,10 +84,8 @@ function nlsolve!(u::AbstractVector, f::AbstractSemidiscreteFunction, cache::New
         linear_solver_cache.isfresh = true # Notify linear solver that we touched the system matrix
 
         residualnorm = residual_norm(cache, f)
-        if residualnorm < cache.parameters.tol
-            if cache.iter == 1
-                push!(Θks, 0.0)
-            end
+        if residualnorm < cache.parameters.tol && cache.iter > 0
+            push!(Θks, 0.0)
             break
         elseif cache.iter > cache.parameters.max_iter
             push!(Θks,Inf)
@@ -116,17 +114,20 @@ function nlsolve!(u::AbstractVector, f::AbstractSemidiscreteFunction, cache::New
             else
                 push!(Θks, Θk)
             end
+            # Try to prevent oversolving when we really just wanted to force the solve to happen once.
+            if cache.iter == 1 && residualnormprev < eps(T) && residualnorm < eps(T) && incrementnorm < eps(T) && incrementnormprev < eps(T)
+                break
+            end
             if cache.parameters.enforce_monotonic_convergence && Θk ≥ 1.0
                 @debug "Newton-Raphson diverged. Aborting. ||r|| = $residualnorm" _group=:nlsolve
                 return false
             end
-
-            # Late out on second iteration
-            if residualnorm < cache.parameters.tol
-                break
-            end
         end
 
+        # if incrementnorm < cache.parameters.tol
+        #     break
+        # end
+
         residualnormprev  = residualnorm
         incrementnormprev = incrementnorm
     end