From f826906e89adfb650141d067c51651178245f144 Mon Sep 17 00:00:00 2001 From: Christian Merdon Date: Wed, 16 Aug 2023 13:04:48 +0200 Subject: [PATCH] small improvements --- examples/Example235_StokesIteratedPenalty.jl | 15 ++----- examples/Example240_SVRTEnrichment.jl | 4 +- examples/Example280_CompressibleStokes.jl | 42 +++++++++++++------- src/common_operators/linear_operator.jl | 2 +- src/solvers.jl | 27 +++---------- 5 files changed, 40 insertions(+), 50 deletions(-) diff --git a/examples/Example235_StokesIteratedPenalty.jl b/examples/Example235_StokesIteratedPenalty.jl index 94b0fc39..603fa9b1 100644 --- a/examples/Example235_StokesIteratedPenalty.jl +++ b/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 diff --git a/examples/Example240_SVRTEnrichment.jl b/examples/Example240_SVRTEnrichment.jl index 9a9268f5..bcf7f47a 100644 --- a/examples/Example240_SVRTEnrichment.jl +++ b/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) diff --git a/examples/Example280_CompressibleStokes.jl b/examples/Example280_CompressibleStokes.jl index effb8d37..954405b6 100644 --- a/examples/Example280_CompressibleStokes.jl +++ b/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,7 +269,13 @@ 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]) @@ -270,7 +283,8 @@ function prepare_data!(; M = 1, c = 1, μ = 1 ) 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 \ No newline at end of file diff --git a/src/common_operators/linear_operator.jl b/src/common_operators/linear_operator.jl index 4eddca6a..81dd6a7e 100644 --- a/src/common_operators/linear_operator.jl +++ b/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) diff --git a/src/solvers.jl b/src/solvers.jl index e880a0f2..70d6f648 100644 --- a/src/solvers.jl +++ b/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