small improvements

examples/Example235_StokesIteratedPenalty.jl

@@ -94,18 +94,9 @@ function main(; Plotter = nothing, λ = 1e4, μ = 1.0, nrefs = 5, kwargs...)
     FES = [FESpace{FETypes[1]}(xgrid), FESpace{FETypes[2]}(xgrid)]
     sol = FEVector([FES[1],FES[2]]; tags = [u,p])
 
-    oldsol = deepcopy(sol)
-    sol, SCu = solve(PDu, [FES[1]]; return_config = true, init = sol, kwargs...)
-    sol, SCp = solve(PDp, [FES[2]]; return_config = true, init = sol, is_linear = true, kwargs...)
-    difference = norm(oldsol.entries - sol.entries)
-    nits = 1
-    while difference > 1e-8
-        nits += 1
-        copyto!(oldsol.entries, sol.entries)
-        sol, SCu = solve(PDu, [FES[1]], SCu; return_config = true, init = sol, kwargs...)
-        sol, SCp = solve(PDp, [FES[2]], SCp; return_config = true, init = sol, is_linear = true, kwargs...)
-        difference = norm(oldsol.entries - sol.entries)
-    end
+    SC1 = SolverConfiguration(PDu, [FES[1]]; init = sol, maxiterations = 1, target_residual = 1e-8, constant_matrix = true, kwargs...)
+    SC2 = SolverConfiguration(PDp, [FES[2]]; init = sol, maxiterations = 1, target_residual = 1e-8, constant_matrix = true, kwargs...)
+    sol, nits = iterate_until_stationarity([SC1, SC2]; init = sol, kwargs...)
     @info "converged after $nits iterations"
 
     ## error calculation

examples/Example240_SVRTEnrichment.jl

@@ -147,11 +147,11 @@ function main(; nrefs = 4, order = 2, Plotter = nothing, enrich = true, reduce =
                uR => enrich ? FES_enrich : nothing)
 
     ## define operators
-    assign_operator!(PD, LinearOperator(rhs!, [id(u)]; bonus_quadorder = 5, kwargs...)) 
+    assign_operator!(PD, LinearOperator(rhs!, [id(u)]; bonus_quadorder = bonus_quadorder, kwargs...)) 
     assign_operator!(PD, BilinearOperator(kernel_stokes_standard!, [grad(u), id(p)]; kwargs...))  
     if enrich
         if reduce
-            @time if order == 1
+            if order == 1
                 @info "... preparing condensation of RT0 dofs"
                 AR = FEMatrix(FES_enrich)
                 BR = FEMatrix(FES[p], FES_enrich)

examples/Example280_CompressibleStokes.jl

@@ -57,10 +57,10 @@ using Symbolics
 ## everything is wrapped in a main function
 ## testcase = 1 : well-balanced test (stratified no-flow over mountain)
 ## testcase = 2 : vortex example (ϱu is div-free p7 vortex)
-function main(; testcase = 1, nrefs = 4, M = 1, c = 1, target_residual = 1e-10, Plotter = nothing, reconstruct = true, μ = 1, order = 1, kwargs...)
+function main(; testcase = 1, nrefs = 4, M = 1, c = 1, ufac = 100, laplacian_in_rhs = false, maxsteps = 5000, target_residual = 1e-10, Plotter = nothing, reconstruct = true, μ = 1, order = 1, kwargs...)
 
 	## load data for testcase
-    grid_builder, kernel_gravity!, u!, ∇u!, ϱ! = load_testcase_data(testcase; nrefs = nrefs, M = M, c = c, μ = μ)
+    grid_builder, kernel_gravity!, kernel_rhs!, u!, ∇u!, ϱ!, τfac = load_testcase_data(testcase; laplacian_in_rhs = laplacian_in_rhs, nrefs = nrefs, M = M, c = c, μ = μ, ufac = ufac)
     xgrid = grid_builder(nrefs)
 
     ## define unknowns
@@ -83,10 +83,13 @@ function main(; testcase = 1, nrefs = 4, M = 1, c = 1, target_residual = 1e-10,
     assign_operator!(PD, BilinearOperator([grad(u)]; factor = μ, store = true, kwargs...))
     assign_operator!(PD, LinearOperator([div(u)], [id(ϱ)]; factor = c, kwargs...))
     assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:4, kwargs...))
+    if kernel_rhs! !== nothing
+        assign_operator!(PD, LinearOperator(kernel_rhs!, [id_u]; factor = 1, store = true, bonus_quadorder = 2*order, kwargs...))
+    end
     assign_operator!(PD, LinearOperator(kernel_gravity!, [id_u], [id(ϱ)]; factor = 1, bonus_quadorder = 3*order, kwargs...))
 
     ## FVM for continuity equation
-	τ = order * μ / (order*M*c) # time step for pseudo timestepping
+	τ = order * μ / (order*M*sqrt(τfac)) # time step for pseudo timestepping
     PDT = ProblemDescription("continuity equation")
     assign_unknown!(PDT, ϱ)    
     if order > 1
@@ -122,9 +125,9 @@ function main(; testcase = 1, nrefs = 4, M = 1, c = 1, target_residual = 1e-10,
         NDofs[lvl] = length(sol.entries)
 
         ## solve the two problems iteratively [1] >> [2] >> [1] >> [2] ...
-        SC1 = SolverConfiguration(PD, FES; init = sol, maxiterations = 1, target_residual = target_residual, constant_matrix = true, kwargs...)
-        SC2 = SolverConfiguration(PDT, FES; init = sol, maxiterations = 1, target_residual = target_residual, kwargs...)
-        sol, nits = iterate_until_stationarity([SC1, SC2]; init = sol, kwargs...)
+        SC1 = SolverConfiguration(PD, [FES[1]]; init = sol, maxiterations = 1, target_residual = target_residual, constant_matrix = true, kwargs...)
+        SC2 = SolverConfiguration(PDT, [FES[2]]; init = sol, maxiterations = 1, target_residual = target_residual, kwargs...)
+        sol, nits = iterate_until_stationarity([SC1, SC2]; maxsteps = maxsteps, init = sol, kwargs...)
 
         ## caculate error
         error = evaluate(ErrorIntegratorExact, sol)
@@ -188,7 +191,7 @@ function standard_gravity!(result, ϱ, qpinfo)
     result[2] = -ϱ[1]
 end
 
-function load_testcase_data(testcase::Int = 1; nrefs = 1, M = 1, c = 1, μ = 1)
+function load_testcase_data(testcase::Int = 1; laplacian_in_rhs = true, nrefs = 1, M = 1, c = 1, μ = 1, ufac = 100)
     if testcase == 1
         grid_builder = (nref) -> simplexgrid(Triangulate;
                     points = [0 0; 0.2 0; 0.3 0.2; 0.45 0.05; 0.55 0.35; 0.65 0.2; 0.7 0.3; 0.8 0; 1 0; 1 1 ; 0 1]',
@@ -203,7 +206,7 @@ function load_testcase_data(testcase::Int = 1; nrefs = 1, M = 1, c = 1, μ = 1)
         M_exact = integrate(xgrid, ON_CELLS, (result, qpinfo) -> (result[1] = exp(-qpinfo.x[2]/c);), 1; quadorder = 20)
         area = sum(xgrid[CellVolumes])
         ϱ1!(result, qpinfo) = (result[1] = exp(-qpinfo.x[2]/c)/(M_exact/area);)
-        return grid_builder, standard_gravity!, u1!, ∇u1!, ϱ1!
+        return grid_builder, standard_gravity!, nothing, u1!, ∇u1!, ϱ1!, 1
     elseif testcase == 2
         grid_builder = (nref) -> simplexgrid(Triangulate;
                     points = [0 0; 1 0; 1 1 ; 0 1]',
@@ -215,7 +218,7 @@ function load_testcase_data(testcase::Int = 1; nrefs = 1, M = 1, c = 1, μ = 1)
 
         xgrid = grid_builder(3)
         M_exact = integrate(xgrid, ON_CELLS, (result, qpinfo) -> (result[1] = exp(-qpinfo.x[1]^3/(3*c));), 1; quadorder = 20)
-        ϱ_eval, g_eval, u_eval, ∇u_eval = prepare_data!(; M = M_exact, c = c, μ = μ)
+        ϱ_eval, g_eval, f_eval, u_eval, ∇u_eval = prepare_data!(; laplacian_in_rhs = laplacian_in_rhs, M = M_exact, c = c, μ = μ, ufac = ufac)
         ϱ2!(result, qpinfo) = (result[1] = ϱ_eval(qpinfo.x[1], qpinfo.x[2]);)
 
         M_exact = integrate(xgrid, ON_CELLS, ϱ2!, 1)
@@ -225,14 +228,18 @@ function load_testcase_data(testcase::Int = 1; nrefs = 1, M = 1, c = 1, μ = 1)
             result .= input[1] * g_eval(qpinfo.x[1], qpinfo.x[2])
         end
 
+        function kernel_rhs!(result, qpinfo)
+            result .= f_eval(qpinfo.x[1], qpinfo.x[2])
+        end
+
         u2!(result, qpinfo) = (result .= u_eval(qpinfo.x[1], qpinfo.x[2]);)
         ∇u2!(result, qpinfo) = (result .= ∇u_eval(qpinfo.x[1], qpinfo.x[2]);)
-        return grid_builder, kernel_gravity!, u2!, ∇u2!, ϱ2!
+        return grid_builder, kernel_gravity!, f_eval === nothing ? nothing : kernel_rhs!, u2!, ∇u2!, ϱ2!, ufac
     end
 end
 
 ## exact data for testcase 2 computed by Symbolics
-function prepare_data!(; M = 1, c = 1, μ = 1 )
+function prepare_data!(; M = 1, c = 1, μ = 1, ufac = 100, laplacian_in_rhs = true)
 
 	@variables x y
 	dx = Differential(x)
@@ -243,7 +250,7 @@ function prepare_data!(; M = 1, c = 1, μ = 1 )
 
 	## stream function ξ
 	## sucht that ϱu = curl ξ
-	ξ = x^2*y^2*(x-1)^2*(y-1)^2 
+	ξ = x^2*y^2*(x-1)^2*(y-1)^2 * ufac
 
 	∇ξ = Symbolics.gradient(ξ, [x,y])
 
@@ -262,15 +269,22 @@ function prepare_data!(; M = 1, c = 1, μ = 1 )
 
 	## gravity ϱg = - Δu + ϱ∇log(ϱ)
 
-	g = - μ*Δu/ϱ + Symbolics.gradient(log(ϱ), [x,y]) 
+    if laplacian_in_rhs
+        f = - μ*Δu
+        g = Symbolics.gradient(log(ϱ), [x,y]) 
+	else
+        g = - μ*Δu/ϱ + Symbolics.gradient(log(ϱ), [x,y]) 
+        f = 0
+    end
 
 	#Δu = Symbolics.derivative(∇u[1,1], [x]) + Symbolics.derivative(∇u[2,2], [y])
 
 	ϱ_eval = build_function(ϱ, x, y, expression = Val{false})
 	u_eval = build_function(u, x, y, expression = Val{false})
 	∇u_eval = build_function(∇u_reshaped, x, y, expression = Val{false})
 	g_eval = build_function(g, x, y, expression = Val{false})
+	f_eval = build_function(f, x, y, expression = Val{false})
 
-    return ϱ_eval, g_eval[1], u_eval[1], ∇u_eval[1]
+    return ϱ_eval, g_eval[1], f == 0 ? nothing : f_eval[1], u_eval[1], ∇u_eval[1]
 end
 end

src/common_operators/linear_operator.jl

@@ -79,7 +79,7 @@ function LinearOperator(kernel, u_test, ops_test, u_args, ops_args; Tv = Float64
     return LinearOperator{Tv, typeof(u_test[1]), typeof(kernel), typeof(storage)}(u_test, ops_test, u_args, ops_args, kernel, [[zeros(Tv, 0, 0, 0)]], [[zeros(Tv, 0, 0, 0)]], nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, storage, parameters)
 end
 
-function LinearOperator(kernel::Function, u_test::Array{<:Unknown,1}, ops_test::Array{DataType,1}; Tv = Float64, kwargs...)
+function LinearOperator(kernel::Function, u_test, ops_test::Array{DataType,1}; Tv = Float64, kwargs...)
     parameters=Dict{Symbol,Any}( k => v[1] for (k,v) in default_linop_kwargs())
     _update_params!(parameters, kwargs)
     @assert length(u_test) == length(ops_test)

src/solvers.jl

@@ -305,7 +305,7 @@ function CommonSolve.solve(PD::ProblemDescription, FES::Array, SC = nothing; unk
 end
 
 
-function iterate_until_stationarity(SCs::Array{<:SolverConfiguration,1}, FES = nothing; maxsteps = 1000, tolerance_change = 1e-10, init = nothing, unknowns = [SC.PD.unknowns for SC in SCs], kwargs... )
+function iterate_until_stationarity(SCs::Array{<:SolverConfiguration,1}, FES = nothing; maxsteps = 1000, init = nothing, unknowns = [SC.PD.unknowns for SC in SCs], kwargs... )
 
     PDs::Array{ProblemDescription,1} = [SC.PD for SC in SCs]
     nPDs = length(PDs)
@@ -331,22 +331,6 @@ function iterate_until_stationarity(SCs::Array{<:SolverConfiguration,1}, FES = n
         sol = FEVector(FES; tags = all_unknowns)
     end
 
-
-    ## init solver configurations
-    time = @elapsed begin
-        allocs = @allocated begin
-            for p = 1 : nPDs
-                SCs[p] = SolverConfiguration(PDs[p], unknowns[p], FES[p]; init = sol, kwargs...)
-                if SCs[p].parameters[:verbosity] > 0
-                    @info ".... init solver configurations \n"
-                end
-                if SCs[p].parameters[:show_config]
-                    @info "\n$(SCs[p])"
-                end
-            end
-        end
-    end
-
     @info "SOLVING iteratively $([PD.name for PD in PDs])
             unknowns = $([[uj.name for uj in u] for u in unknowns])"
       #      fetypes = $(["$(get_FEType(FES[j]))" for j = 1 : length(unknowns)])
@@ -385,9 +369,8 @@ function iterate_until_stationarity(SCs::Array{<:SolverConfiguration,1}, FES = n
 
     alloc_factor = 1024^2
 
-    @printf " init time | allocs = %.2f s | %.2f MiB\n" time allocs/alloc_factor
-    time_final = time
-    allocs_final = allocs
+    time_final = 0
+    allocs_final = 0
     nlres = 1.1e30
     linres = 1.1e30
     converged = zeros(Bool, nPDs)
@@ -479,7 +462,9 @@ function iterate_until_stationarity(SCs::Array{<:SolverConfiguration,1}, FES = n
 
                 time_solve = @elapsed begin
                     allocs_solve = @allocated begin
-                        linsolve.A = A.entries.cscmatrix
+                        if !SC.parameters[:constant_matrix] || it == 1
+                            linsolve.A = A.entries.cscmatrix
+                        end
                         linsolve.b = b.entries
 
                         ## solve