residual in solve until stationarity solver are saved now (#115)

* residual in solve until stationarity solver are saved now and can be extracted from statistics, printing now verbosity dependent

* use residuals in tests of Example280

Author: chmerdon

CHANGELOG.md

@@ -1,6 +1,12 @@
 # CHANGES
 
 
+## v.1.10.5
+
+### Fixed
+  - linear and nonlinear residuals in solve_until_stationarity are saved in SolverConfiguration.statistics
+
+
 ## v.1.10.4
 
 ### Fixed

Project.toml

@@ -1,6 +1,6 @@
 name = "ExtendableFEM"
 uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed"
-version = "1.10.4"
+version = "1.10.5"
 authors = ["Christian Merdon <merdon@wias-berlin.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de>"]
 
 [deps]

examples/Example280_CompressibleStokes.jl

@@ -130,7 +130,7 @@ function main(;
     EnergyIntegrator = ItemIntegrator(energy_kernel!, [id(u)]; resultdim = 1, quadorder = 2 * (order + 1), kwargs...)
     ErrorIntegratorExact = ItemIntegrator(exact_error!(u!, ∇u!, ϱ!), [id(u), grad(u), id(ϱ)]; resultdim = 9, quadorder = 2 * (order + 1), kwargs...)
     NDofs = zeros(Int, nrefs)
-    Results = zeros(Float64, nrefs, 5)
+    Results = zeros(Float64, nrefs, 7)
 
     sol = nothing
     xgrid = nothing
@@ -152,13 +152,20 @@ function main(;
         SC2 = SolverConfiguration(PDT; init = sol, maxiterations = 1, target_residual = target_residual, kwargs...)
         sol, nits = iterate_until_stationarity([SC1, SC2]; energy_integrator = EnergyIntegrator, maxsteps = maxsteps, init = sol, kwargs...)
 
+        residual_momentum = residual(SC1)
+        residual_continuity = residual(SC2)
+        @info "final residual momentum = $(residual_momentum)"
+        @info "final residual continuity = $(residual_continuity)"
+
         ## calculate error
         error = evaluate(ErrorIntegratorExact, sol)
         Results[lvl, 1] = sqrt(sum(view(error, 1, :)) + sum(view(error, 2, :)))
         Results[lvl, 2] = sqrt(sum(view(error, 3, :)) + sum(view(error, 4, :)) + sum(view(error, 5, :)) + sum(view(error, 6, :)))
         Results[lvl, 3] = sqrt(sum(view(error, 7, :)))
         Results[lvl, 4] = sqrt(sum(view(error, 8, :)) + sum(view(error, 9, :)))
         Results[lvl, 5] = nits
+        Results[lvl, 6] = residual_momentum
+        Results[lvl, 7] = residual_continuity
 
         ## print results
         print_convergencehistory(NDofs[1:lvl], Results[1:lvl, :]; X_to_h = X -> X .^ (-1 / 2), ylabels = ["|| u - u_h ||", "|| ∇(u - u_h) ||", "|| ϱ - ϱ_h ||", "|| ϱu - ϱu_h ||", "#its"], xlabel = "ndof")
@@ -328,6 +335,8 @@ end
 generateplots = ExtendableFEM.default_generateplots(Example280_CompressibleStokes, "example280.png") #hide
 function runtests() #hide
     Results, plt = main(; nrefs = 2) #hide
+    @test Results[end, 6] <= 1.0e-11 #hide
+    @test Results[end, 7] <= 1.0e-11 #hide
     @test Results[end, 1] ≈ 6.732891488265023e-7 #hide
     return nothing #hide
 end #hide

src/solvers.jl

@@ -868,11 +868,14 @@ function iterate_until_stationarity(
     allocs_final = 0
     nlres = 1.1e30
     linres = 1.1e30
+    verbosity = maximum([SC.parameters[:verbosity] for SC in SCs])
     converged = zeros(Bool, nPDs)
     it::Int = 0
     while (it < maxsteps) && (any(converged .== false))
         it += 1
-        @printf "%5d\t" it
+        if verbosity > -2
+            @printf "%5d\t" it
+        end
         copyto!(init.entries, sol.entries)
         allocs_assembly = 0
         time_assembly = 0
@@ -884,6 +887,7 @@ function iterate_until_stationarity(
             A = As[p]
             PD = PDs[p]
             SC = SCs[p]
+            stats = SC.statistics
             residual = residuals[p]
             maxits = SC.parameters[:maxiterations]
             nltol = SC.parameters[:target_residual]
@@ -923,7 +927,10 @@ function iterate_until_stationarity(
                         end
                     end
                     nlres = norm(residual.entries)
-                    @printf "\tres[%d] = %.2e" p nlres
+                    if verbosity > -2
+                        @printf "\tres[%d] = %.2e" p nlres
+                    end
+                    push!(stats[:nonlinear_residuals], nlres)
                 end
                 time_final += time_assembly + time_solve_init
                 allocs_final += allocs_assembly + allocs_solve_init
@@ -944,15 +951,18 @@ function iterate_until_stationarity(
                 allocs_final += allocs_solve
                 time_solve += time_solve_init
                 allocs_solve += allocs_solve_init
-                if SC.parameters[:verbosity] > -1
+                if verbosity > -1
                     @printf " (%.3e)" linres
                 end
+                push!(stats[:linear_residuals], linres)
             end # nonlinear iterations subproblem
         end
 
         if energy_integrator !== nothing
             error = evaluate(energy_integrator, sol)
-            @printf "   energy = %.3e" sum([sum(view(error, j, :)) for j in 1:size(error, 1)])
+            if verbosity > -1
+                @printf "   energy = %.3e" sum([sum(view(error, j, :)) for j in 1:size(error, 1)])
+            end
         end
         @printf "\n"
     end