From ec4f98a0546d986456926885c6ee7d18294a87b9 Mon Sep 17 00:00:00 2001 From: Alexander Korzec Date: Thu, 25 Jun 2026 20:47:51 +0200 Subject: [PATCH 1/3] Added .gitignore, Complemented Project.toml, Minor code corrections --- .gitignore | 7 +++++++ Project.toml | 11 +++++++++++ scripts/poisson.jl | 4 ++-- src/modelproblems/logpoisson_dual.jl | 2 +- src/modelproblems/solvers_logpoisson_primal.jl | 1 - src/modelproblems/solvers_poisson_primal.jl | 1 - 6 files changed, 21 insertions(+), 5 deletions(-) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..d6a81a5 --- /dev/null +++ b/.gitignore @@ -0,0 +1,7 @@ +# Folders +.CondaPkg +data +plots + +# Files +Manifest.toml \ No newline at end of file diff --git a/Project.toml b/Project.toml index c9e43d0..aa75ac4 100644 --- a/Project.toml +++ b/Project.toml @@ -4,9 +4,12 @@ version = "1.0.1" authors = ["Christian Merdon ", "Martin Eigel "] [deps] +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +DrWatson = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1" ExtendableFEM = "a722555e-65e0-4074-a036-ca7ce79a4aed" ExtendableFEMBase = "12fb9182-3d4c-4424-8fd1-727a0899810c" ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8" @@ -14,9 +17,12 @@ ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3" GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" +LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +PythonPlot = "274fc56d-3b97-40fa-a1cd-1b4a50311bf9" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" @@ -24,10 +30,13 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [compat] Aqua = "0.8" +CairoMakie = "0.15.9" +DataFrames = "1.8.1" Distributed = "1.9" Distributions = "0.25.120" DocStringExtensions = "0.9.4" DoubleFloats = "1.4.3" +DrWatson = "2.19.1" ExplicitImports = "1" ExtendableFEM = "1.10" ExtendableFEMBase = "1.5.1" @@ -36,8 +45,10 @@ ExtendableSparse = "2.0.1" GridVisualize = "1.12" IterativeSolvers = "0.9.2" Krylov = "0.10.1" +LaTeXStrings = "1.4.0" LinearAlgebra = "1.9" Printf = "1.9" +PythonPlot = "1.0.6" Random = "1.9" SparseArrays = "1.9" SpecialFunctions = "2.5.1" diff --git a/scripts/poisson.jl b/scripts/poisson.jl index 1b4bfb3..ed50f18 100644 --- a/scripts/poisson.jl +++ b/scripts/poisson.jl @@ -3,7 +3,7 @@ ([source code](SOURCE_URL)) runs AFEM loop for stochastic Poisson problem - + usage: - run experiment: run(; problem = problem, kwargs...) - load results : show_results(; kwargs...) @@ -484,7 +484,7 @@ function produce_plots(; for o in 1:length(order), d in 1:length(decay) basedata["order"] = order[o] basedata["decay"] = decay[d] - data[o, d], ~ = produce_or_load(main, basedata, filename = filename, force = force) + data[o, d], ~ = produce_or_load(_main, basedata, filename = filename, force = force) xgrid = data[o, d]["solution"].FES_space[1].xgrid repair_grid!(xgrid) xgrid[BFaceRegions] .= 1 diff --git a/src/modelproblems/logpoisson_dual.jl b/src/modelproblems/logpoisson_dual.jl index c44e052..a7f12cf 100644 --- a/src/modelproblems/logpoisson_dual.jl +++ b/src/modelproblems/logpoisson_dual.jl @@ -39,7 +39,7 @@ function deterministic_problem(::Type{LogTransformedPoissonProblemDual}, C::Abst end PD = ProblemDescription("log-transformed Poisson problem (dual)") - p = Unknown("p"; name = "stress p = -∇ũ - ∇aũ = - epx(-a) ∇u") + p = Unknown("p"; name = "stress p = -∇ũ - ∇aũ = - exp(-a) ∇u") u = Unknown("u"; name = "potential ũ = exp(-a) u") assign_unknown!(PD, p) assign_unknown!(PD, u) diff --git a/src/modelproblems/solvers_logpoisson_primal.jl b/src/modelproblems/solvers_logpoisson_primal.jl index c78d794..0258ef4 100644 --- a/src/modelproblems/solvers_logpoisson_primal.jl +++ b/src/modelproblems/solvers_logpoisson_primal.jl @@ -177,7 +177,6 @@ function solve_logpoisson_primal_full!(SolutionSGFEM::SGFEVector, A, N0, N, b, G M::Int = length(N) # size(G,1) / nmodes FES = SolutionSGFEM.FES_space[1] - nmodes = num_multiindices(SolutionSGFEM) bigFES = [FES for j in 1:nmodes] x::Vector{Float64} = zeros(Float64, 2) diff --git a/src/modelproblems/solvers_poisson_primal.jl b/src/modelproblems/solvers_poisson_primal.jl index 543d496..de0677c 100644 --- a/src/modelproblems/solvers_poisson_primal.jl +++ b/src/modelproblems/solvers_poisson_primal.jl @@ -176,7 +176,6 @@ function solve_full_primal!(SolutionSGFEM::SGFEVector, A0, A, b, G, nmodes, rhsf FES = SolutionSGFEM.FES_space[1] @show FES multi_indices = SolutionSGFEM.TB.multi_indices - nmodes = num_multiindices(SolutionSGFEM) bigFES = [FES for j in 1:nmodes] x::Vector{Float64} = zeros(Float64, 2) From 791a25ed0e62a54509c7da370b045937e9a6ce47 Mon Sep 17 00:00:00 2001 From: Alexander Korzec Date: Thu, 25 Jun 2026 20:57:51 +0200 Subject: [PATCH 2/3] Created PoC version of StokesSimple and Stokes scripts, Defined StokesProblemPrimal, Implemented full matrix and matrix-free solvers for StokesProblemPrimal, Implemented pressure error computation with more generic calculate_sampling_error alternative. Unfortunately, the gradient does not converge for the rhs in the Stokes script --- scripts/stokes.jl | 584 +++++++++++++++++++++ scripts/stokes_simple.jl | 86 +++ src/ExtendableASGFEM.jl | 13 +- src/estimate.jl | 29 + src/modelproblems/modelproblems.jl | 1 + src/modelproblems/solvers_stokes_primal.jl | 266 ++++++++++ src/modelproblems/stokes_primal.jl | 93 ++++ src/sampling_error.jl | 172 +++++- 8 files changed, 1238 insertions(+), 6 deletions(-) create mode 100644 scripts/stokes.jl create mode 100644 scripts/stokes_simple.jl create mode 100644 src/modelproblems/solvers_stokes_primal.jl create mode 100644 src/modelproblems/stokes_primal.jl diff --git a/scripts/stokes.jl b/scripts/stokes.jl new file mode 100644 index 0000000..bf7e111 --- /dev/null +++ b/scripts/stokes.jl @@ -0,0 +1,584 @@ +#= + +([source code](SOURCE_URL)) + +runs AFEM loop for stochastic Stokes problem + +usage: +- run experiment: run(; problem = problem, kwargs...) +- load results : show_results(; kwargs...) +- produce plots : produce_plots(; kwargs...) + +possible values for problem are +- StokesProblemPrimal = Stokes problem with linear coefficient ν + +=# + +module Stokes + +using CairoMakie +using DataFrames +using DelimitedFiles +using DrWatson +using ExtendableASGFEM +using ExtendableFEM +using ExtendableFEMBase +using ExtendableGrids +using GridVisualize +using LaTeXStrings +using Pkg +using UnicodePlots + +function f!(result, qpinfo) + x = qpinfo.x[1] + y = qpinfo.x[2] + result[1] = 5(x^4) + 12(x^2) * ((-1 + x)^2) * (-1 + y) + 12(x^2) * ((-1 + x)^2) * y + 4(x^2) * y * ((-1 + y)^2) + 4(x^2) * (-1 + y) * (y^2) + 4((-1 + x)^2) * (-1 + y) * (y^2) + 4((-1 + x)^2) * y * ((-1 + y)^2) + 16x * (-1 + x) * (-1 + y) * (y^2) + 16x * (-1 + x) * y * ((-1 + y)^2) + result[2] = 5(y^4) - 4x * ((-1 + x)^2) * ((-1 + y)^2) - 4(-1 + x) * (x^2) * ((-1 + y)^2) - 4x * ((-1 + x)^2) * (y^2) - 4(-1 + x) * (x^2) * (y^2) - 16(-1 + x) * (x^2) * y * (-1 + y) - 16x * ((-1 + x)^2) * y * (-1 + y) - 12x * (y^2) * ((-1 + y)^2) - 12(-1 + x) * (y^2) * ((-1 + y)^2) + #result[1] = 1000 * (1 - x) * x * y * (1 - y) + #result[2] = 1000 * (1 - x) * x * y * (1 - y) + #result[1] = sin(x) * sin(y) * cos(pi * x / 2) * cos(pi * y / 2) + #ξ = x^2 * (x - 1)^2 * y^2 * (y - 1)^2 + #result[1] = - (1 / 6) * x^3 + #result[2] = (1 / 6) * y^3 + return nothing +end + +function filename(data; folder = "data", add = "", makepath = false) + problem = data["problem"] + decay = data["decay"] + mean = data["mean"] + order = data["order"] + domain = data["domain"] + θ_stochastic = data["θ_stochastic"] + θ_spatial = data["θ_spatial"] + tail_extension = data["tail_extension"] + maxdofs = data["maxdofs"] + filename = "$(folder)/$(problem)/$(domain)/" + if makepath + mkpath(filename) + end + filename *= add * "order=$(order)_maxdofs=$(maxdofs)_decay=$(decay)_mean=$(mean)_θ=($(θ_spatial),$(θ_stochastic))_tail=$(tail_extension)" + return filename +end + +default_args = Dict( + "initial_refs" => 1, + "domain" => "square", + "order" => 1, + "problem" => StokesProblemPrimal, + "C" => StochasticCoefficientCosinus, + "decay" => 1000, + "mean" => 1.0, + "maxm" => 150, + "bonus_quadorder_a" => 2, + "f" => f!, + "bonus_quadorder_f" => 0, + "θ_stochastic" => 0.5, + "θ_spatial" => 0.5, + "factor_tail" => 1, + "tail_extension" => [10, 2], # for [0] mode, and all others + "maxdofs" => 1.0e4, + "initial_modes" => [[0]], + "nsamples" => 150, + "use_iterative_solver" => true, + "Plotter" => nothing, +) + +function run(; force = false, kwargs...) + data = deepcopy(default_args) + + for (k, v) in kwargs + data[String(k)] = v + end + Plotter = data["Plotter"] + + data, ~ = produce_or_load(_main, data, filename = filename, force = force) + + results = data["results"] + sol = data["solution"] + multi_indices = data["multi_indices"] + repair_grid!(sol.FES_space[1].xgrid) + if Plotter === nothing + for j in 1:num_multiindices(sol) + @show extrema(view(sol[j])) + println(stdout, unicode_scalarplot(sol[j]; title = "MI $(multi_indices[j])")) + end + else + plot_modes(sol; Plotter = Plotter, ncols = 4) + end + + @info "final multi_indices = " + for m in multi_indices + println("$m") + end + @show results + return data +end + +function show_results(; force = false, mode_history_up_to_level = 15, kwargs...) + data = deepcopy(default_args) + + for (k, v) in kwargs + data[String(k)] = v + end + + data, ~ = produce_or_load(_main, data, filename = filename, force = force) + + results = data["results"] + nmodes = results[!, "nmodes"] + nlevels = length(nmodes) + + MI = data["multi_indices"] + for j in 1:min(mode_history_up_to_level, nlevels) + istart = j == 1 ? 1 : nmodes[j - 1] + 1 + iend = nmodes[j] + if iend >= istart + @info "NEW MODES on level $j : " + maxlength = maximum([findlast(!=(0), MI[k]) for k in istart:iend]) + if maxlength === nothing + maxlength = 1 + end + for k in istart:iend + @info MI[k][1:maxlength] + end + end + end + + return data +end + +function _main( + data = nothing; + debug = false, + plot_solution = false, + Plotter = nothing, + kwargs... + ) + if isnothing(data) + data = deepcopy(default_args) + end + + for (k, v) in kwargs + data[String(k)] = v + end + + problem = data["problem"] + nrefs = data["initial_refs"] + decay = data["decay"] + mean = data["mean"] + τ = 0.9 + maxm = data["maxm"] + C = data["C"](; τ = τ, decay = decay, mean = mean, maxm = maxm) + θ_stochastic = data["θ_stochastic"] + θ_spatial = data["θ_spatial"] + multi_indices = Array{Array{Int, 1}, 1}(data["initial_modes"]) + nsamples = data["nsamples"] + f! = data["f"] + bonus_quadorder_a = data["bonus_quadorder_a"] + bonus_quadorder_f = data["bonus_quadorder_f"] + factor_tail = data["factor_tail"] + maxdofs = data["maxdofs"] + order = data["order"] + use_iterative_solver = data["use_iterative_solver"] + + prepare_multi_indices!(multi_indices) + @info "initial multi_indices = $(multi_indices)" + + dim = 2 + domain = data["domain"] + if dim == 1 + xgrid = uniform_refine(reference_domain(Edge1D), nrefs) + elseif dim == 2 + if domain == "square" + xgrid = uniform_refine(grid_unitsquare(Triangle2D), nrefs) + elseif domain == "lshape" + xgrid = uniform_refine(grid_lshape(Triangle2D), nrefs) + else + @error "domain $domain not known" + end + else + @error "dimension $dim not available yet" + end + + df = DataFrame( + ndofs_space = zeros(Int, 0), + nmodes = zeros(Int, 0), + exact_error_stress = zeros(Float64, 0), + exact_error_u = zeros(Float64, 0), + exact_error_p = zeros(Float64, 0), + exact_error_stress_uni = zeros(Float64, 0), + exact_error_u_uni = zeros(Float64, 0), + exact_error_p_uni = zeros(Float64, 0), + estimate = zeros(Float64, 0), + estimate_active = zeros(Float64, 0), + estimate_tail = zeros(Float64, 0), + time_solve = zeros(Float64, 0), + time_estimate = zeros(Float64, 0), + ) + + sol = nothing + lvl = 0 + + while true + push!(df, [0 0 0 0 0 0 0 0 0 0 0 0 0]) + + lvl += 1 + @info "LEVEL $lvl" + prepare_multi_indices!(multi_indices) + M = maximum(length.(multi_indices)) + OBType = LegendrePolynomials + ansatz_deg = maximum([maximum(multi_indices[k]) for k in 1:length(multi_indices)]) + 4 + TensorBasis = TensorizedBasis(OBType, M, ansatz_deg, 2 * ansatz_deg, 2 * ansatz_deg, multi_indices = multi_indices) + if debug + @info C, TensorBasis + end + + if num_nodes(xgrid) < 200 + println(stdout, unicode_gridplot(xgrid)) + end + + FETypes = (H1BR{2}, L2P0{1}) + FES = [FESpace{FETypes[1]}(xgrid), FESpace{FETypes[2]}(xgrid)] + unames = ["u", "p"] + + sol = SGFEVector(FES, TensorBasis; active_modes = 1:length(multi_indices), unames = unames) + + time_solve = @elapsed bdofs = solve!(problem, sol, C; rhs = (f!), bonus_quadorder_a = bonus_quadorder_a, bonus_quadorder_f = bonus_quadorder_f, use_iterative_solver = use_iterative_solver, debug = debug) + df[lvl, :time_solve] = time_solve + + if isnothing(Plotter) && plot_solution + for j in 1:num_multiindices(sol) + if debug + @show extrema(view(sol[j])) + end + println(stdout, unicode_scalarplot(sol[j]; title = "MI $(multi_indices[j])")) + end + elseif plot_solution + plot_modes(sol; Plotter = Plotter, ncols = 4) + end + weightederrorH1, weightederrorL2u, weightederrorL2p, uniformerrorH1, uniformerrorL2u, uniformerrorL2p = calculate_sampling_error_2( + sol, C; problem = problem, metrics_configurations = stokes_metrics_configuration, rhs = (f!), order = order + 1, nsamples, debug + ) + + tail_extension = data["tail_extension"] + time_estimate = @elapsed η4modes, η4cells, multi_indices_extended = estimate(problem, sol, C; rhs = f!, bonus_quadorder = max(bonus_quadorder_a, bonus_quadorder_f), tail_extension = tail_extension) + for (multi_index, error) in zip(multi_indices_extended, η4modes) + @info "mode = $(multi_index) | error = $(error)" + end + df[lvl, :time_estimate] = time_estimate + + inactive_else, inactive_bnd, inactive_bnd2, active_bnd, active_int = classify_modes(multi_indices_extended, multi_indices_extended[1:TensorBasis.nmodes]) + + sorted_id = sortperm(η4modes; rev = true) + EstimatorInterior = 0 + EstimatorBoundary = 0 + actives = union(active_int, active_bnd) + for j in 1:length(multi_indices_extended) + print("$(multi_indices_extended[sorted_id[j]]) | ") + if sorted_id[j] in inactive_else + print(" [inactive far]") + elseif sorted_id[j] in inactive_bnd + print("[inactive bnd1]") + elseif sorted_id[j] in inactive_bnd2 + print("[inactive bnd2]") + elseif sorted_id[j] in active_int + print(" [active int] ") + elseif sorted_id[j] in active_bnd + print(" [active bnd] ") + end + println(" | $(η4modes[sorted_id[j]])") + if sorted_id[j] in actives + EstimatorInterior += η4modes[sorted_id[j]]^2 + else + EstimatorBoundary += η4modes[sorted_id[j]]^2 + end + end + @show actives + + @info "FINAL RESULTS + || u - u_h || (w,u) = $(sqrt.(weightederrorL2u)), $(sqrt.(uniformerrorL2u)) + || ∇(u - u_h) || (w,u) = $(sqrt.(weightederrorH1)), $(sqrt.(uniformerrorH1)) + || p - p_h || (w,u) = $(sqrt.(weightederrorL2p)), $(sqrt.(uniformerrorL2p)) + estimator_total = $(sqrt(sum(η4modes .^ 2))) + estimator_space = $(sqrt(EstimatorInterior)) + estimator_stoch = $(sqrt(EstimatorBoundary))" + + df[lvl, :ndofs_space] = sum([sol.FES_space[j].ndofs for j in 1:length(sol.FES_space)]) + df[lvl, :nmodes] = sol.TB.nmodes + df[lvl, :exact_error_stress] = sqrt(weightederrorH1[end]) + df[lvl, :exact_error_u] = sqrt(weightederrorL2u[end]) + df[lvl, :exact_error_p] = sqrt(weightederrorL2p[end]) + df[lvl, :exact_error_stress_uni] = sqrt(uniformerrorH1[end]) + df[lvl, :exact_error_u_uni] = sqrt(uniformerrorL2u[end]) + df[lvl, :exact_error_p_uni] = sqrt(uniformerrorL2p[end]) + df[lvl, :estimate] = sqrt(sum(η4modes .^ 2)) + df[lvl, :estimate_active] = sqrt(EstimatorInterior) + df[lvl, :estimate_tail] = sqrt(EstimatorBoundary) + data["results"] = df + data["solution"] = sol + data["multi_indices"] = multi_indices + + if length(sol.entries) >= maxdofs + break + end + + if EstimatorInterior > factor_tail * EstimatorBoundary + println("...spatial refinement") + if θ_spatial >= 1 + xgrid = uniform_refine(xgrid) + else + facemarker = bulk_mark(xgrid, view(sum(view(η4cells, :, actives), dims = 2), :), θ_spatial; indicator_AT = ON_CELLS) + xgrid = RGB_refine(xgrid, facemarker) + end + else + println("...stochastic refinement") + accumulated_stochastic_error = 0 + for j in 1:length(multi_indices_extended) + if (sorted_id[j] in inactive_bnd || sorted_id[j] in inactive_bnd2) + percentage = η4modes[sorted_id[j]]^2 / EstimatorBoundary + println("... adding mode $(multi_indices_extended[sorted_id[j]]) with error $(η4modes[sorted_id[j]]) (=$(Float16(percentage * 100))% of tail error)") + accumulated_stochastic_error += η4modes[sorted_id[j]]^2 + push!(multi_indices, multi_indices_extended[sorted_id[j]]) + if accumulated_stochastic_error >= θ_stochastic * EstimatorBoundary + break + end + end + end + end + + @info "MODE-CLASSIFICATION + active_interior = $(active_int) + active_boundary = $(active_bnd) + inactive_bnd = $(inactive_bnd) + inactive_bnd2 = $(inactive_bnd2) + inactive_else = $(inactive_else)" + end + + filename_params = filename(data; add = "", folder = "data", makepath = true) * ".txt" + + data["version"] = pkgversion(ExtendableASGFEM) + writedlm(filename_params, data, "=") + + return data +end + +function produce_plots(; + order = default_args["order"], + decay = default_args["decay"], + legend_position = :lb, + xscale = log10, + yscale = log10, + force = false, + scaling_factor = true, + markersize = 10, + plotM = 80, + maxdegree = 12, + template = :convergence, + cols = nothing, + xlabel = "ndofs", + show_optimal_rate = false, + kwargs... + ) + colors = Makie.wong_colors() + linestyles = [:solid, :dot, :dashdot, :dashdotdot] + + basedata = deepcopy(default_args) + for (k, v) in kwargs + basedata[String(k)] = v + end + + if template == :convergence + cols = isnothing(cols) ? [:estimate, :exact_error_stress, :exact_error_p] : cols + elseif template == :convergence_with_L2 + cols = isnothing(cols) ? [:estimate, :exact_error_stress, :exact_error, :exact_error_p] : cols + elseif template == :active_vs_tail + cols = isnothing(cols) ? [:estimate_active, :estimate_tail] : cols + elseif template == :dofs + cols = isnothing(cols) ? [:ndofs_space, :nmodes, :ndofs_all] : cols + xscale = nothing + legend_position = :lt + xlabel = "level" + else + @error "template = $template not implemented" + end + + if typeof(order) <: Real + order = [order] + end + if typeof(decay) <: Real + decay = [decay] + end + data = Array{Dict{String, Any}, 2}(undef, length(order), length(decay)) + for o in 1:length(order), d in 1:length(decay) + basedata["order"] = order[o] + basedata["decay"] = decay[d] + data[o, d], ~ = produce_or_load(_main, basedata, filename = filename, force = force) + xgrid = data[o, d]["solution"].FES_space[1].xgrid + repair_grid!(xgrid) + xgrid[BFaceRegions] .= 1 + + gplt = GridVisualize.gridplot(xgrid; Plotter = CairoMakie, linewidth = 1, colorbar = :none, legend = :none) + filename_plot = filename(basedata; add = "grid_", folder = "plots", makepath = true) * ".png" + CairoMakie.save(filename_plot, gplt) + @info "grid plot for $cols saved under $filename_plot" + end + + for d in 1:length(decay) + ylabel = "" + if length(cols) == 1 + if cols[1] == :exact_error_u + ylabel = L"|| u - u_h ||_{L^2}" + elseif cols[1] == :exact_error_stress + ylabel = L"|| u - u_h ||_A" + elseif col[1] == :exect_error_p + ylabel = L"|| p - p_h ||_{L^2}" + elseif cols[1] == :estimate + ylabel = L"\eta" + elseif cols[1] == :estimate_tail + ylabel = L"\eta_\text{tail}" + elseif cols[1] == :estimate_active + ylabel = L"\eta_\text{active}" + end + end + title = "" + f_conv = Figure(fontsize = 18, size = (900, 600)) + if xscale !== nothing + ax_conv = Axis( + f_conv[1, 1], title = title, xlabel = xlabel, ylabel = ylabel, yscale = yscale, xscale = xscale, + xminorticksvisible = true, xminorgridvisible = true, yminorgridvisible = true, yminorticksvisible = true, + xminorticks = IntervalsBetween(10), yminorticks = IntervalsBetween(10) + ) + else + ax_conv = Axis( + f_conv[1, 1], title = title, xlabel = xlabel, ylabel = ylabel, yscale = yscale, + xminorticksvisible = true, xminorgridvisible = true, yminorgridvisible = true, yminorticksvisible = true, + xminorticks = IntervalsBetween(10), yminorticks = IntervalsBetween(10) + ) + end + + @info "producing convergence history..." + for o in 1:length(order) + df = data[o, d]["results"] + nlevels = size(df, 1) + ndofs = nothing + ndofs_space = nothing + ndofs_stochastic = nothing + for c in 1:length(cols) + color = colors[o] + ndofs_space = df[!, :ndofs_space] + ndofs_stochastic = df[!, :nmodes] + if template == :dofs + ndofs = 1:length(ndofs_space) + else + ndofs = df[!, :ndofs_space] .* df[!, :nmodes] + end + if cols[c] == :ndofs_all + plotdata = df[!, :ndofs_space] .* df[!, :nmodes] + else + plotdata = df[!, cols[c]] + end + + label = "" + if length(cols) > 1 + if cols[c] == :exact_error_u + label = L"|| u - u_h ||_{L^2}" + elseif cols[c] == :exact_error_stress + label = L"|| u - u_h ||_A" + elseif cols[c] == :estimate + if scaling_factor !== nothing + plotdata .*= scaling_factor + end + label = L"\eta" + elseif cols[c] == :exact_error_p + label = L"|| p - p_h ||_{L^2}" + elseif cols[c] == :estimate_tail + if scaling_factor !== nothing + plotdata .*= scaling_factor + end + label = L"\eta(\partial_h \Lambda)" + elseif cols[c] == :estimate_active + if scaling_factor !== nothing + plotdata .*= scaling_factor + end + label = L"\eta(\Lambda)" + elseif cols[c] == :ndofs_space + label = L"| V_h |" + elseif cols[c] == :nmodes + label = L"| \Lambda |" + elseif cols[c] == :ndofs_all + label = L"| \mathcal{V}_h |" + end + end + label = L"%$label order = %$(order[o])" + scatterlines!(ndofs, plotdata, label = label, color = color, markersize = markersize, linewidth = 2, linestyle = linestyles[c]) + end + if show_optimal_rate + label = L"\mathcal{O}(N_\text{space}^{%$(-order[o]/2)})" + scatterlines!(ndofs, ndofs .^ (-order[o] / 2), label = label, color = :gray, markersize = 0, linewidth = 2, linestyle = linestyles[length(cols) + 1]) + end + end + axislegend(ax_conv, position = legend_position, merge = true, labelsize = 16, orientation = :horizontal, nbanks = length(order) * length(cols) + 1) + + basedata["order"] = order + basedata["decay"] = decay[d] + filename_plot = filename(basedata; add = "$(template)_", folder = "plots", makepath = true) * ".png" + CairoMakie.save(filename_plot, f_conv) + @info "convergence plot for $cols saved under $filename_plot" + end + + ## bar plot + for o in 1:length(order) + max4mode = zeros(length(decay), plotM) + cat = zeros(Int, plotM * length(decay)) + height = zeros(Float64, plotM * length(decay)) + groups = zeros(Int, plotM * length(decay)) + for d in 1:length(decay) + multi_indices = data[o, d]["multi_indices"] + nmodes = length(multi_indices) + M = min(maximum(length.(multi_indices)), plotM) + max4mode[d, 1:M] = [maximum([multi_indices[j][k] for j in 1:nmodes]) for k in 1:M] + cat[d:length(decay):end] .= 1:plotM + groups[d:length(decay):end] .= d + height[d:length(decay):end] .= max4mode[d, :] + end + + + # Plot + replace!(height, 0 => 0.1) + yticks = (1:maxdegree, ["$(Int(j))" for j in 1:maxdegree]) + xticks = (1:plotM, [j % 2 == 1 ? "$(Int(j))" : "" for j in 1:plotM]) + f_bar = Figure(fontsize = 12, size = (900, 250)) + Axis(f_bar[1, 1], limits = (0, plotM, 0, maxdegree), yticks = yticks, xticks = xticks, xlabel = "stochastic dimension", ylabel = "max degree", xticklabelrotation = pi / 2) + CairoMakie.barplot!( + cat, height, + dodge = AbstractVector{Integer}(groups), + color = colors[groups] + ) + + # Legend + labels = ["decay $d" for d in decay] + elements = [PolyElement(polycolor = colors[i]) for i in 1:length(labels)] + title = "" + + Legend(f_bar[1, 2], elements, labels, title) + + basedata["decay"] = decay + basedata["order"] = order[o] + filename_plot = filename(basedata; add = "barplot_", folder = "plots", makepath = true) * ".png" + CairoMakie.save(filename_plot, f_bar) + @info "modes bar plot for $cols saved under $filename_plot" + end + return +end + +function repair_grid!(xgrid::ExtendableGrid) + xgrid.components[CellGeometries] = VectorOfConstants{ElementGeometries, Int}(xgrid.components[CellGeometries][1], num_cells(xgrid)) + xgrid.components[FaceGeometries] = VectorOfConstants{ElementGeometries, Int}(xgrid.components[FaceGeometries][1], length(xgrid.components[FaceGeometries])) + xgrid.components[BFaceGeometries] = VectorOfConstants{ElementGeometries, Int}(xgrid.components[BFaceGeometries][1], length(xgrid.components[BFaceGeometries])) + + xgrid.components[UniqueCellGeometries] = Vector{ElementGeometries}([xgrid.components[CellGeometries][1]]) + xgrid.components[UniqueFaceGeometries] = Vector{ElementGeometries}([xgrid.components[FaceGeometries][1]]) + return xgrid.components[UniqueBFaceGeometries] = Vector{ElementGeometries}([xgrid.components[BFaceGeometries][1]]) +end + +end diff --git a/scripts/stokes_simple.jl b/scripts/stokes_simple.jl new file mode 100644 index 0000000..9257b2a --- /dev/null +++ b/scripts/stokes_simple.jl @@ -0,0 +1,86 @@ +module StokesSimple + +using ExtendableASGFEM +using ExtendableFEM +using ExtendableFEMBase +using ExtendableGrids +using GridVisualize + +function f!(result, qpinfo) + x = qpinfo.x[1] + y = qpinfo.x[2] + result[1] = 5(x^4) + 12(x^2) * ((-1 + x)^2) * (-1 + y) + 12(x^2) * ((-1 + x)^2) * y + 4(x^2) * y * ((-1 + y)^2) + 4(x^2) * (-1 + y) * (y^2) + 4((-1 + x)^2) * (-1 + y) * (y^2) + 4((-1 + x)^2) * y * ((-1 + y)^2) + 16x * (-1 + x) * (-1 + y) * (y^2) + 16x * (-1 + x) * y * ((-1 + y)^2) + result[2] = 5(y^4) - 4x * ((-1 + x)^2) * ((-1 + y)^2) - 4(-1 + x) * (x^2) * ((-1 + y)^2) - 4x * ((-1 + x)^2) * (y^2) - 4(-1 + x) * (x^2) * (y^2) - 16(-1 + x) * (x^2) * y * (-1 + y) - 16x * ((-1 + x)^2) * y * (-1 + y) - 12x * (y^2) * ((-1 + y)^2) - 12(-1 + x) * (y^2) * ((-1 + y)^2) + return nothing + #result[1] = 5(x^4) - 0.001(-12(x^2) * ((-1 + x)^2) * (-1 + y) - 12(x^2) * ((-1 + x)^2) * y - 4(x^2) * y * ((-1 + y)^2) - 4(x^2) * (-1 + y) * (y^2) - 4((-1 + x)^2) * (-1 + y) * (y^2) - 4((-1 + x)^2) * y * ((-1 + y)^2) - 16x * (-1 + x) * (-1 + y) * (y^2) - 16x * (-1 + x) * y * ((-1 + y)^2)) + #result[2] = 5(y^4) - 0.001(4x * ((-1 + x)^2) * ((-1 + y)^2) + 4(-1 + x) * (x^2) * ((-1 + y)^2) + 4x * ((-1 + x)^2) * (y^2) + 4(-1 + x) * (x^2) * (y^2) + 16(-1 + x) * (x^2) * y * (-1 + y) + 16x * ((-1 + x)^2) * y * (-1 + y) + 12x * (y^2) * ((-1 + y)^2) + 12(-1 + x) * (y^2) * ((-1 + y)^2)) + #result[1] = 1000 * (1 - x) * x * y * (1 - y) + #result[2] = 1000 * (1 - x) * x * y * (1 - y) +end + +function main(; + problem = StokesProblemPrimal, + nrefs = 3, # number of uniform refinements of the initial grid + order = 2, # polynomial order of the FEspaces + decay = 2.0, # decay factor for the random coefficient + mean = 1, # mean value of coefficient + domain = "square", # domain, e.g., "square" or "lshape" + initial_modes = [[0], [1, 0], [0, 1], [2, 0], [0, 0, 1]], # initial multi-indices for stochastic basis + (f!) = (f!), # right-hand side function + use_iterative_solver = true, # use iterative solver ? (otherwise direct) + nsamples = 20, + Plotter = nothing, + debug = false, + ) + ## prepare stochastic coefficient + C = StochasticCoefficientCosinus(; τ = 0.9, decay = decay, mean = mean) + + ## prepare grid + xgrid = if domain == "square" + uniform_refine(grid_unitsquare(Triangle2D), nrefs) + elseif domain == "lshape" + uniform_refine(grid_lshape(Triangle2D), nrefs) + else + error("unknown domain: $domain") + end + + ## prepare stochastic basis + multi_indices = Array{Array{Int, 1}, 1}(initial_modes) + prepare_multi_indices!(multi_indices) + M = maximum(length.(multi_indices)) + OBType = LegendrePolynomials + ansatz_deg = maximum([maximum(multi_indices[k]) for k in 1:length(multi_indices)]) + 4 + TensorBasis = TensorizedBasis(OBType, M, ansatz_deg, 2 * ansatz_deg, 2 * ansatz_deg, multi_indices = multi_indices) + + ## prepare FE spaces + FETypes = (H1BR{2}, L2P0{1}) + FES = [FESpace{FETypes[1]}(xgrid), FESpace{FETypes[2]}(xgrid)] + unames = ["u", "p"] + + ## create solution vector + sol = SGFEVector(FES, TensorBasis; active_modes = 1:length(multi_indices), unames = unames) + + ## solve problem + @info "Solving..." + solve!(problem, sol, C; rhs = (f!), use_iterative_solver = use_iterative_solver) + + ## plot solution + if !isnothing(Plotter) + p = plot_modes(sol; Plotter = Plotter, ncols = 4) + display(p) + end + + ## compute exact error (by MC sampling) + weightederrorH1, weightederrorL2u, weightederrorL2p, uniformerrorH1, uniformerrorL2u, uniformerrorL2p = calculate_sampling_error_2( + sol, C; problem = problem, metrics_configurations = stokes_metrics_configuration, rhs = (f!), order = order + 1, nsamples, debug + ) + + @info "RESULTS + || ∇(u-u_h) || (w,u) = $(sqrt(weightederrorH1[end])), $(sqrt(uniformerrorH1[end])) + || u - u_h || (w,u) = $(sqrt(weightederrorL2u[end])), $(sqrt(uniformerrorL2u[end])) + || p - p_h || (w,u) = $(sqrt(weightederrorL2p[end])), $(sqrt(uniformerrorL2p[end]))" + + return sol +end + +end diff --git a/src/ExtendableASGFEM.jl b/src/ExtendableASGFEM.jl index cbb9e41..5719f2d 100644 --- a/src/ExtendableASGFEM.jl +++ b/src/ExtendableASGFEM.jl @@ -5,13 +5,13 @@ using Distributions: Distributions, Normal, Uniform, dim, pdf using DocStringExtensions: DocStringExtensions, TYPEDEF, TYPEDSIGNATURES using ExtendableFEM: ExtendableFEM, BilinearOperator, FaceInterpolator, HomogeneousBoundaryData, ItemIntegrator, ItemIntegratorDG, - L2NormIntegrator, LinearOperator, ProblemDescription, - Unknown, assemble!, assign_operator!, assign_unknown!, - div, grad, id, jump, plot + L2NormIntegrator, LinearOperator, ProblemDescription, Reconstruct, Unknown, + ZeroMeanValueRestriction, apply, assemble!, assign_operator!, assign_restriction!, + assign_unknown!, div, grad, id, jump, plot using ExtendableFEMBase: ExtendableFEMBase, BFaceDofs, CellDofs, Divergence, FEEvaluator, FEMatrix, FESpace, FEVector, - FEVectorBlock, Gradient, H1P1, H1Pk, HDIVRTk, - Laplacian, QuadratureRule, _addnz, addblock!, + FEVectorBlock, Gradient, H1P1, H1Pk, H1P2B, L2P1, HDIVRT0, HDIVRT1, HDIVRTk, + H1BR, Identity, L2P0, Laplacian, QuadratureRule, _addnz, addblock!, addblock_matmul!, eval_febe!, fill!, get_ndofs, get_polynomialorder, unicode_scalarplot, update_basis! @@ -85,9 +85,12 @@ export deterministic_problem, deterministic_problem2 export LogTransformedPoissonProblemPrimal export LogTransformedPoissonProblemDual export PoissonProblemPrimal +export StokesProblemPrimal include("sampling_error.jl") export calculate_sampling_error +export calculate_sampling_error_2 +export u_with_stress_metric_configuration, stokes_metrics_configuration include("estimate.jl") export estimate diff --git a/src/estimate.jl b/src/estimate.jl index aadadb3..48197d4 100644 --- a/src/estimate.jl +++ b/src/estimate.jl @@ -379,6 +379,35 @@ function estimate(::Type{PoissonProblemPrimal}, sol::SGFEVector, C::AbstractStoc return eta4modes, eta4cell, multi_indices_extended #, f4modes end +# TODO: Implement StokesProblemPrimal a posteriori estimator +# This is a dummy estimator which enforces repeated spatial refinement +function estimate(::Type{StokesProblemPrimal}, sol::SGFEVector, C::AbstractStochasticCoefficient; rhs = nothing, bonus_quadorder = 1, tail_extension = 5) + FES = sol.FES_space[1] + xgrid = FES.xgrid + ncells = num_cells(xgrid) + TB = sol.TB + nmodes = TB.nmodes + multi_indices = TB.multi_indices + + multi_indices_extended = add_boundary_modes(deepcopy(multi_indices); tail_extension = tail_extension) + + nmodes_extended = length(multi_indices_extended) + + inactive_else, inactive_bnd, inactive_bnd2, active_bnd, active_int = classify_modes(multi_indices_extended, multi_indices_extended[1:nmodes]) + + eta4modes = ones(Float64, nmodes_extended) + eta4cell = ones(Float64, ncells, nmodes_extended) + + actives = union(active_int, active_bnd) + for j in 1:length(nmodes_extended) + if j in actives + eta4modes[j] = 10000 + end + end + + return eta4modes, eta4cell, multi_indices_extended +end + # function estimate(::Type{LogTransformedPoissonProblemDual}, sol::SGFEVector, C::AbstractStochasticCoefficient; problem = LogTransformedPoissonProblemPrimal, rhs = nothing, bonus_quadorder = 0) diff --git a/src/modelproblems/modelproblems.jl b/src/modelproblems/modelproblems.jl index bb0afd3..812c8be 100644 --- a/src/modelproblems/modelproblems.jl +++ b/src/modelproblems/modelproblems.jl @@ -29,3 +29,4 @@ end include("logpoisson_primal.jl") include("logpoisson_dual.jl") include("poisson_primal.jl") +include("stokes_primal.jl") diff --git a/src/modelproblems/solvers_stokes_primal.jl b/src/modelproblems/solvers_stokes_primal.jl new file mode 100644 index 0000000..1fd648f --- /dev/null +++ b/src/modelproblems/solvers_stokes_primal.jl @@ -0,0 +1,266 @@ +#################### +### SOLVER STUFF ### +#################### + +struct StokesPrimal{Tv, MT, VT, GT} + A0::MT + A::Vector{MT} + B::MT + G::GT + bdofs::Vector{Int} + nmodes::Int + vsize::Array{Int, 1} +end + +struct StokesPrimalPreconditioner{Tv, FAC} + LUS::FAC + DA::Array{Tv, 1} + temp::Array{Tv, 1} + bdofs::Vector{Int} + nmodes::Int + vsize::Array{Int, 1} +end + +function stokesPrimalPreconditioner(A0::ExtendableSparseMatrix{Tv, Ti}, B::ExtendableSparseMatrix{Tv, Ti}, bdofs, nmodes, vsize) where {Tv, Ti} + #for dof in bdofs + # A0[dof, dof] = 1.0e60 + #end + #flush!(A0) + + DA::Array{Tv, 1} = zeros(Tv, size(A0, 1)) + for j in 1:length(DA) + DA[j] = A0[j, j] + end + + # compute S = B inv(A_diag) B' + S = ExtendableSparseMatrix{Tv, Ti}(size(B, 2), size(B, 2)) + cscmat::SparseMatrixCSC{Tv, Ti} = B.cscmatrix + rows::Array{Ti, 1} = rowvals(cscmat) + valsB::Array{Tv, 1} = cscmat.nzval + value::Tv = 0 + row::Ti = 0 + for i in 1:size(B, 2), j in 1:size(B, 2) + for r in nzrange(cscmat, i), r2 in nzrange(cscmat, j) + if rows[r] == rows[r2] + row = rows[r] + value = valsB[r] * valsB[r2] / DA[row] + _addnz(S, i, j, value, 1) + end + end + end + + # compute LU factorisation of S + flush!(S) + LUS = lu(S.cscmatrix) + #LUA = lu(A0.cscmatrix) + + # temporary storage array for solver + temp = zeros(Tv, size(B, 2)) + + return StokesPrimalPreconditioner{Tv, typeof(LUS)}(LUS, DA, temp, bdofs, nmodes, vsize) +end + +@inline LinearAlgebra.ldiv!(C::StokesPrimalPreconditioner, b) = ldiv!(b, C, b) +@inline function LinearAlgebra.ldiv!(y, C::StokesPrimalPreconditioner{Tv, FAC}, b) where {Tv, FAC} + a::Int = 0 + c::Int = 0 + DA::Array{Tv, 1} = C.DA + temp::Array{Tv, 1} = C.temp + nmodes::Int = C.nmodes + vsize::Array{Int, 1} = C.vsize + for mu in 1:nmodes + # upper left block of preconditioner (I ⊗ A_diag) + a = (mu - 1) * vsize[1] + 1 + c = mu * vsize[1] + for i in a:c + y[i] = b[i] / DA[i - a + 1] + end + a = nmodes * vsize[1] + (mu - 1) * vsize[2] + 1 + c = nmodes * vsize[1] + mu * vsize[2] + if y !== b + ldiv!(view(y, a:c), C.LUS, view(b, a:c)) + else + ldiv!(temp, C.LUS, view(b, a:c)) + y[a:c] .= temp + end + end + + return y +end + +@inline function LinearAlgebra.:\(C::StokesPrimalPreconditioner, b) + y = zero(b) + ldiv!(y, C, b) + return y +end + +function LinearAlgebra.mul!(Ax::Vector{Tv}, S::StokesPrimal{Tv, MT, VT, GT}, x) where {Tv, MT, VT, GT} + fill!(Ax, 0) + g::Tv = 0 + vsize::Array{Int, 1} = S.vsize + nmodes::Int = S.nmodes + bdofs::Vector{Int} = S.bdofs + G::GT = S.G + A::Vector{MT} = S.A + M::Int = length(A) # size(G,1) / nmodes + a::Int = 0 + b::Int = 0 + a2::Int = 0 + b2::Int = 0 + A0::MT = S.A0 + B::MT = S.B + for mu in 1:nmodes + # deterministic part + a = (mu - 1) * vsize[1] + 1 + b = mu * vsize[1] + a2 = a + b2 = b + addblock_matmul!(view(Ax, a:b), A0[1, 1], view(x, a2:b2)) + a2 = nmodes * vsize[1] + (mu - 1) * vsize[2] + 1 + b2 = nmodes * vsize[1] + mu * vsize[2] + addblock_matmul!(view(Ax, a:b), B[1, 1], view(x, a2:b2)) + a2 = a + b2 = b + a = nmodes * vsize[1] + (mu - 1) * vsize[2] + 1 + b = nmodes * vsize[1] + mu * vsize[2] + addblock_matmul!(view(Ax, a:b), B[1, 1], view(x, a2:b2); transposed = true) + + # stochastic part + a = (mu - 1) * vsize[1] + 1 + b = mu * vsize[1] + for nu in 1:nmodes, e in 1:M + g = G[(e - 1) * nmodes + mu, nu] + if abs(g) > 1.0e-12 + a2 = (nu - 1) * vsize[1] + 1 + b2 = nu * vsize[1] + addblock_matmul!(view(Ax, a:b), A[e][1, 1], view(x, a2:b2); factor = g) + end + end + + for dof in bdofs + Ax[a + dof - 1] = 0 + end + end + + return nothing +end + +Base.eltype(S::StokesPrimal) = typeof(S).parameters[1] +Base.size(S::StokesPrimal) = S.nmodes .* (size(S.A0.entries) .+ size(S.B.entries)[2]) + +function solve_stokes_primal!(SolutionSGFEM::SGFEVector, A0, A, B, b0, G, nmodes, bfac; atol = 1.0e-14, rtol = 1.0e-14) + ## create fullmatrix-free matrix evaluator + @info "Solving StochasticFEM iteratively and matrix-free (ndofs = $(length(SolutionSGFEM.entries)))..." + vsize = [SolutionSGFEM[1].FES.ndofs, SolutionSGFEM[nmodes + 1].FES.ndofs] + ## boundary data + bfacedofs = SolutionSGFEM.FES_space[1][BFaceDofs] + nbfaces = num_sources(bfacedofs) + bdofs = [] + for bface in 1:nbfaces + append!(bdofs, view(bfacedofs, :, bface)) + end + bdofs = unique(bdofs) + + S = StokesPrimal{eltype(G), typeof(A0), typeof(b0), typeof(G)}(A0, A, B, G, bdofs, nmodes, vsize) + @info "...initializing Preconditioner" + @time P = stokesPrimalPreconditioner(A0.entries, B.entries, bdofs, nmodes, vsize) + + ## right-hand side + b = deepcopy(SolutionSGFEM) + fill!(b.entries, 0) + addblock!(b[1], b0[1]; factor = bfac) + for m in 1:nmodes + for dof in bdofs + b[m][dof] = 0 + end + end + + ## solve + @info "...starting right-conditioned GMRES" + x, history = Krylov.gmres(S, b.entries, SolutionSGFEM.entries; ldiv = true, atol = atol, rtol = rtol, M = P) + SolutionSGFEM.entries .= x + @show history + + ## Pressure mean von pressure modes abziehen, damit Vorfaktor wie bei use_iterative_solver=false ist. + ## https://wias-pdelib.github.io/ExtendableFEM.jl/stable/module_examples/Example252_NSEPlanarLatticeFlow/ + xgrid = SolutionSGFEM[1].FES.xgrid + for i in 1:nmodes + pintegrate = ItemIntegrator([id(1)]) + pmean = sum(ExtendableFEM.evaluate(pintegrate, [SolutionSGFEM.FEVectorBlocks[nmodes + i]])) / sum(xgrid[CellVolumes]) + view(SolutionSGFEM.FEVectorBlocks[nmodes + i]) .-= pmean + end + + ## check residual + Ax = zero(SolutionSGFEM.entries) + mul!(Ax, S, SolutionSGFEM.entries) + @info "solver residual = $(sqrt(sum((Ax - b.entries) .^ 2)))" + + return bdofs +end + +function solve_stokes_primal_full!(SolutionSGFEM::SGFEVector, A0, A, B, b0, G, nmodes, rhsfac) + M::Int = length(A) # size(G,1) / nmodes + FES = SolutionSGFEM.FES_space + multi_indices = SolutionSGFEM.TB.multi_indices + nmodes = num_multiindices(SolutionSGFEM) + bigFES = Array{FESpace{Float64, Int32}, 1}([FES[1] for j in 1:nmodes]) + append!(bigFES, [FES[2] for j in 1:nmodes]) + + bigS = FEMatrix(bigFES) + bigb = FEVector(bigFES) + + for j in 1:nmodes + addblock!(bigS[j, j], A0[1, 1]) + addblock!(bigS[j, nmodes + j], B[1, 1]) + addblock!(bigS[nmodes + j, j], B[1, 1]; transpose = true) + end + + ## right-hand side + addblock!(bigb[1], b0[1]) + + g::Float64 = 0 + for j in 1:nmodes, k in 1:nmodes + for e in 1:M + g = G[(e - 1) * nmodes + j, k] # ⟨ ξ_m ψ_mi(j) ψ_mi(k) ⟩ + if abs(g) > 1.0e-12 + addblock!(bigS[j, k], A[e][1, 1]; factor = g) + end + end + end + flush!(bigS.entries) + + ## boundary data + bfacedofs = FES[1][BFaceDofs] + nbfaces = num_sources(bfacedofs) + bdofs = [] + for bface in 1:nbfaces + append!(bdofs, view(bfacedofs, :, bface)) + end + unique!(bdofs) + + for m in 1:nmodes + for dof in bdofs + bigS[m, m][dof, dof] = 1.0e60 + bigb[m][dof] = 0 + end + end + flush!(bigS.entries) + + @info "Solving StochasticFEM with full matrix..." + SolutionSGFEM.entries .= bigS.entries \ bigb.entries + + xgrid = SolutionSGFEM[1].FES.xgrid + for i in 1:nmodes + pintegrate = ItemIntegrator([id(1)]) + pmean = sum(ExtendableFEM.evaluate(pintegrate, [SolutionSGFEM.FEVectorBlocks[nmodes + i]])) / sum(xgrid[CellVolumes]) + view(SolutionSGFEM.FEVectorBlocks[nmodes + i]) .-= pmean + end + + residual = bigS.entries * SolutionSGFEM.entries .- bigb.entries + for m in 1:nmodes + residual[FES[1].ndofs * (m - 1) .+ bdofs] .= 0 + end + println("linear residual = $(sqrt(sum(residual .^ 2)))") + + return bdofs +end diff --git a/src/modelproblems/stokes_primal.jl b/src/modelproblems/stokes_primal.jl new file mode 100644 index 0000000..1e635ac --- /dev/null +++ b/src/modelproblems/stokes_primal.jl @@ -0,0 +1,93 @@ +""" + +Stokes problem with linear stochastic coefficient ``ν`` that seeks ``(u, p)`` such that + +``-\\nu\\Delta u + \\nabla p = f \\quad \\text{for } (y,x) \\in \\Gamma \\times D`` +``-\\nu\\mathrm{div} u = 0 \\quad \\text{for } (y,x) \\in \\Gamma \\times D`` + +""" +abstract type StokesProblemPrimal <: AbstractModelProblem end + +include("solvers_stokes_primal.jl") + +## deterministic problem description +function deterministic_problem(::Type{StokesProblemPrimal}, C::AbstractStochasticCoefficient, sample_pointer; (get_a!) = (get_a!), rhs = nothing, bonus_quadorder_a = 2, bonus_quadorder_f = 0) + get_ν! = get_a!(C) + + function stokes_kernel!(result, input, qpinfo) + ν = zeros(1) + get_ν!(ν, qpinfo.x, sample_pointer) + result[1] = ν[1] * input[1] - input[5] + result[2] = ν[1] * input[2] + result[3] = ν[1] * input[3] + result[4] = ν[1] * input[4] - input[5] + result[5] = -(input[1] + input[4]) + return nothing + end + + PD = ProblemDescription("Stokes problem (primal)") + u = Unknown("u", name = "velocity", dim = 2) + p = Unknown("p", name = "pressure", dim = 1) + assign_unknown!(PD, u) + assign_unknown!(PD, p) + assign_operator!(PD, BilinearOperator(stokes_kernel!, [grad(u), id(p)]; bonus_quadorder = bonus_quadorder_a)) + if rhs !== nothing + # TODO: Define reconstruct::Bool switch for reconstruction operator + assign_operator!(PD, LinearOperator(rhs, [apply(u, Reconstruct{HDIVRT1{2}, Identity})]; bonus_quadorder = bonus_quadorder_f)) + end + assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:4)) + assign_restriction!(PD, ZeroMeanValueRestriction(p)) + + return PD, u +end + +## solver for stochastic Galerkin problem +function solve!( + ::Type{StokesProblemPrimal}, + sol::SGFEVector, ## target SGFEM vector + C::AbstractStochasticCoefficient; ## stochastic coefficient a + bonus_quadorder_a = 2, ## additional quadrature order for grad(am) + bonus_quadorder_f = 0, ## additional quadrature order for rhs f + rhs = nothing, + debug = false, + use_iterative_solver = true + ) + FES = sol.FES_space + TB = sol.TB + + ## Laplacian: A = (a∇u,∇v) with building blocks Am = (a_m∇u,∇v) + A0 = FEMatrix(FES[1]) + assemble!(A0, BilinearOperator(get_am_x(0, C), [grad(1)], [grad(1)]; bonus_quadorder = bonus_quadorder_a)) + A = [] + + for m in 1:maxlength_multiindices(TB) + Am = FEMatrix(FES[1]) + assemble!(Am, BilinearOperator(get_am_x(m, C), [grad(1)], [grad(1)]; bonus_quadorder = bonus_quadorder_a)) + push!(A, Am) + end + + ## Deterministic B = (∇ ⋅ v, q) term + B = FEMatrix(FES[1], FES[2]) + assemble!(B, BilinearOperator([div(1)], [id(1)]; factor = -1)) # Checken! + + ## assemble right-hand side (f, v) + if rhs !== nothing + b0 = FEVector(FES[1]) + assemble!(b0, LinearOperator(rhs, [apply(1, Reconstruct{HDIVRT0{2}, Identity})]; bonus_quadorder = bonus_quadorder_f)) + else + @error "need right-hand side" + end + + ## get G matrix + G = TB.G + nmodes = num_multiindices(TB) + + ## solve + if use_iterative_solver + @time bdofs = solve_stokes_primal!(sol, A0, A, B, b0, G, nmodes, 1) + else + @time bdofs = solve_stokes_primal_full!(sol, A0, A, B, b0, G, nmodes, 1) + end + + return bdofs +end diff --git a/src/sampling_error.jl b/src/sampling_error.jl index b3510da..444d6b4 100644 --- a/src/sampling_error.jl +++ b/src/sampling_error.jl @@ -1,5 +1,5 @@ ## kernel for L2 error of stress, i.e || grad(u) - grad(u_h) || -function data_error_stress(::Union{Type{<:PoissonProblemPrimal}, Type{<:LogTransformedPoissonProblemPrimal}}, dim, C::AbstractStochasticCoefficient, sample_pointer) +function data_error_stress(::Type{<:PoissonProblemPrimal}, dim, C::AbstractStochasticCoefficient, sample_pointer) function closure(result, input, qpinfo) result[1] = 1 if dim == 1 @@ -27,6 +27,20 @@ function data_error_stress(::Type{<:LogTransformedPoissonProblemDual}, dim, C::A return closure, [id(1), id(2)], [(1, 1), (2, 1)] end +function data_error_stress(::Union{Type{<:StokesProblemPrimal}, Type{<:LogTransformedPoissonProblemPrimal}}, dim, C::AbstractStochasticCoefficient, sample_pointer) + function closure(result, input, qpinfo) + result[1] = 1 + if dim == 1 + result[2] = (input[1] - input[2])^2 + elseif dim == 2 + result[2] = (input[1] - input[3])^2 + (input[2] - input[4])^2 + end + result[1] *= result[2] + return nothing + end + return closure, [grad(1), grad(2)], [(1, 1), (2, 1)] +end + ## kernel for L2 error || u - u_h || function data_error_u(::Union{Type{<:PoissonProblemPrimal}, Type{<:LogTransformedPoissonProblemPrimal}}, dim, C::AbstractStochasticCoefficient, sample_pointer) function closure(result, input, qpinfo) @@ -35,7 +49,24 @@ function data_error_u(::Union{Type{<:PoissonProblemPrimal}, Type{<:LogTransforme end return closure, [id(1), id(2)], [(1, 1), (2, 1)] end + function data_error_u(::Type{<:LogTransformedPoissonProblemDual}, dim, C::AbstractStochasticCoefficient, sample_pointer) + function closure(result, input, qpinfo) + result[1] = (input[1] - input[2])^2 + return nothing + end + return closure, [id(1), id(2)], [(1, 2), (2, 1)] +end + +function data_error_u(::Type{<:StokesProblemPrimal}, dim, C::AbstractStochasticCoefficient, sample_pointer) + function closure(result, input, qpinfo) + result[1] = (input[1] - input[2])^2 + return nothing + end + return closure, [id(1), id(2)], [(1, 1), (2, 1)] +end + +function data_error_p(::Type{<:StokesProblemPrimal}, dim, C::AbstractStochasticCoefficient, sample_pointer) function closure(result, input, qpinfo) result[1] = (input[1] - input[2])^2 return nothing @@ -46,7 +77,33 @@ end FES4sampling(::Type{PoissonProblemPrimal}, dim, xgrid, order) = [FESpace{H1Pk{1, dim, order}}(xgrid)] FES4sampling(::Type{LogTransformedPoissonProblemPrimal}, dim, xgrid, order) = [FESpace{H1Pk{1, dim, order}}(xgrid)] FES4sampling(::Type{LogTransformedPoissonProblemDual}, dim, xgrid, order) = [FESpace{HDIVRTk{dim, order}}(xgrid), FESpace{H1Pk{1, dim, order}}(xgrid; broken = true)] +FES4sampling(::Type{StokesProblemPrimal}, _, xgrid, order) = [FESpace{H1P2B{2, 2}}(xgrid), FESpace{L2P1{1}}(xgrid)] +u_with_stress_metric_configuration = [ + Dict( + "name" => "L2stress", + "closure" => data_error_stress, + ), + Dict( + "name" => "L2u", + "closure" => data_error_u, + ), +] + +stokes_metrics_configuration = [ + Dict( + "name" => "L2stress", + "closure" => data_error_stress, + ), + Dict( + "name" => "L2u", + "closure" => data_error_u, + ), + Dict( + "name" => "L2p", + "closure" => data_error_p, + ), +] """ $(TYPEDSIGNATURES) @@ -217,3 +274,116 @@ function calculate_sampling_error( return totalerrorL2stress_weighted, totalerrorL2u_weighted, totalerrorL2stress_uniform, totalerrorL2u_uniform end + +function calculate_sampling_error_2( + SolutionSGFEM::SGFEVector, + C::AbstractStochasticCoefficient; + metrics_configurations = u_with_stress_metric_configuration, + problem = LogTransformedPoissonProblemPrimal, + bonus_quadorder_a = 10, + bonus_quadorder_f = 0, + order = 2, + rhs = nothing, + dim = size(SolutionSGFEM.FES_space[1].xgrid[Coordinates], 1), + Msamples = maxm(C), + parallel_sampling = true, + dimensionwise_error = true, + energy_norm = true, + debug = false, + nsamples = 100 + ) + + nthreads = Threads.nthreads() + @info "Estimating exact error by MC sampling (with nthreads = $nthreads)" + FES = SolutionSGFEM.FES_space + xgrid = FES[1].xgrid + sol_sgfem = SolutionSGFEM.FEV + TensorBasis = SolutionSGFEM.TB + multi_indices::Array{Array{Int, 1}, 1} = TensorBasis.multi_indices + M::Int = maxlength_multiindices(TensorBasis) + nmodes = num_multiindices(TensorBasis) + + ## generate samples + Samples, weights = sample_distribution(TensorBasis, nsamples; M = Msamples, Mweights = Msamples) + + ## prepare array with deterministic solutions + sol_det = Array{FEVector, 1}(undef, nsamples) + + ## compute deterministic solutions (in parallel) + Threads.@threads for s in 1:nsamples + + ## deterministic problem description + PD, u = deterministic_problem(problem, C, Samples[:, s]; rhs = rhs, bonus_quadorder_a = bonus_quadorder_a, bonus_quadorder_f = bonus_quadorder_f) + + ## solve problem for the current sample + FESSampling = FES4sampling(problem, dim, xgrid, order) + sol_det[s] = ExtendableFEM.solve(PD, FESSampling; verbosity = debug ? 0 : -1, timeroutputs = :none) + print(".") + end + println(" MC samples solved") + + ## compute errors for each sample (sequentially) + M0 = dimensionwise_error ? 0 : M + error_per_run = Dict() + metrics = [] + csample = Samples[:, 1] + metric_names = [configuration["name"] for configuration in metrics_configurations] + for configuration in metrics_configurations + kernel!, input, ids = configuration["closure"](problem, dim, C, csample) + push!( + metrics, + Dict( + "name" => configuration["name"], + "ids" => ids, + "integrator" => ItemIntegrator(kernel!, input; resultdim = dim, quadorder = 2 * order), + "solution" => Array{FEVectorBlock, 1}(undef, 2), + ), + ) + error_per_run[configuration["name"]] = zeros(Float64, M + 1, nsamples) + end + + for s in 1:nsamples + csample = Samples[:, s] + for metric in metrics + for entry in metric["ids"] + if entry[1] == 1 + metric["solution"][1] = sol_det[s][entry[2]] + else + metric["solution"][2] = sol_sgfem[entry[2]] + end + end + end + + for m in M0:M + set_sample!(SolutionSGFEM, view(csample, 1:m)) + for metric in metrics + sol = ExtendableFEM.evaluate(metric["integrator"], metric["solution"]) + error_per_run[metric["name"]][m + 1, s] = sum(view(sol, 1, :)) # Hier Abweichung: Keine Abhängigkeit von energy_norm! + end + end + + @info "SAMPLE $s of $nsamples + sample = $(view(csample, 1:min(Msamples, 2 * M))) + weight = $(weights[s]) + $(join(map(name -> "error($name) = $(error_per_run[name][end, s])", metric_names), "\n "))" + end + + error = Dict() + weightsum = sum(weights) + for name in metric_names + name_weighted = "$(name)_weighted" + name_uniform = "$(name)_uniform" + error[name_weighted] = zeros(Float64, M + 1) + error[name_uniform] = zeros(Float64, M + 1) + for s in 1:nsamples + error[name_uniform] .+= view(error_per_run[name], :, s) + error[name_weighted] .+= view(error_per_run[name], :, s) * weights[s] + end + for m in 1:(M + 1) + error[name_uniform][m] /= nsamples + error[name_weighted][m] /= weightsum + end + end + + return [error["$(name)_weighted"] for name in metric_names]..., [error["$(name)_uniform"] for name in metric_names]... +end From b48506341f515272d1b0cdad330cfb6b4c30fa70 Mon Sep 17 00:00:00 2001 From: Alexander Korzec Date: Thu, 2 Jul 2026 15:45:53 +0200 Subject: [PATCH 3/3] Fixed error calculation for Stokes example and added error calculation for arbitrary dimensions, implemented random constant coefficients, extended Stokes solvers to include boundary conditions --- scripts/stokes.jl | 58 +++++---- scripts/stokes_simple.jl | 68 ++++++----- src/ExtendableASGFEM.jl | 8 +- src/coefficients/coefficients.jl | 2 +- src/coefficients/constant_coefficients.jl | 46 +++++++ src/estimate.jl | 4 +- src/modelproblems/solvers_stokes_primal.jl | 99 +++++++-------- src/modelproblems/stokes_primal.jl | 91 ++++++++++---- src/plots.jl | 1 + src/sampling_error.jl | 133 +++++++++------------ 10 files changed, 305 insertions(+), 205 deletions(-) create mode 100644 src/coefficients/constant_coefficients.jl diff --git a/scripts/stokes.jl b/scripts/stokes.jl index bf7e111..6d9c3b2 100644 --- a/scripts/stokes.jl +++ b/scripts/stokes.jl @@ -1,3 +1,5 @@ +# Stochastisches Beispiel mit Null Randdaten +# Qu & Xu Bsp. 1 nachcoden -> coefficients.jl dazu definieren, an inhomogene Randdaten denken -> Muss implementiert werden #= ([source code](SOURCE_URL)) @@ -27,20 +29,32 @@ using ExtendableGrids using GridVisualize using LaTeXStrings using Pkg +using Symbolics using UnicodePlots -function f!(result, qpinfo) - x = qpinfo.x[1] - y = qpinfo.x[2] - result[1] = 5(x^4) + 12(x^2) * ((-1 + x)^2) * (-1 + y) + 12(x^2) * ((-1 + x)^2) * y + 4(x^2) * y * ((-1 + y)^2) + 4(x^2) * (-1 + y) * (y^2) + 4((-1 + x)^2) * (-1 + y) * (y^2) + 4((-1 + x)^2) * y * ((-1 + y)^2) + 16x * (-1 + x) * (-1 + y) * (y^2) + 16x * (-1 + x) * y * ((-1 + y)^2) - result[2] = 5(y^4) - 4x * ((-1 + x)^2) * ((-1 + y)^2) - 4(-1 + x) * (x^2) * ((-1 + y)^2) - 4x * ((-1 + x)^2) * (y^2) - 4(-1 + x) * (x^2) * (y^2) - 16(-1 + x) * (x^2) * y * (-1 + y) - 16x * ((-1 + x)^2) * y * (-1 + y) - 12x * (y^2) * ((-1 + y)^2) - 12(-1 + x) * (y^2) * ((-1 + y)^2) - #result[1] = 1000 * (1 - x) * x * y * (1 - y) - #result[2] = 1000 * (1 - x) * x * y * (1 - y) - #result[1] = sin(x) * sin(y) * cos(pi * x / 2) * cos(pi * y / 2) - #ξ = x^2 * (x - 1)^2 * y^2 * (y - 1)^2 - #result[1] = - (1 / 6) * x^3 - #result[2] = (1 / 6) * y^3 - return nothing +function prepare_data(C::AbstractStochasticCoefficient, sample_pointer) + @variables ν x y + + ξ = x^2 * (x - 1)^2 * y^2 * (y - 1)^2 + ∇ξ = Symbolics.gradient(ξ, [x, y]) + u = [-∇ξ[2], ∇ξ[1]] + p = x^5 + y^5 - 1 / 3 + ∇u = Symbolics.jacobian(u, [x, y]) + Δu = [ + Symbolics.derivative(∇u[1, 1], x) + Symbolics.derivative(∇u[1, 2], y), + Symbolics.derivative(∇u[2, 1], x) + Symbolics.derivative(∇u[2, 2], y), + ] + ∇p = Symbolics.gradient(p, [x, y]) + f = -ν * Δu + ∇p + eval_f! = build_function(f, ν, x, y, expression = Val{false})[2] + + get_ν! = get_a!(C) + ν = zeros(1) + return function (result, qpinfo) + get_ν!(ν, qpinfo.x, sample_pointer) + eval_f!(result, ν[1], qpinfo.x...) + return nothing + end end function filename(data; folder = "data", add = "", makepath = false) @@ -71,7 +85,7 @@ default_args = Dict( "mean" => 1.0, "maxm" => 150, "bonus_quadorder_a" => 2, - "f" => f!, + "f" => nothing, "bonus_quadorder_f" => 0, "θ_stochastic" => 0.5, "θ_spatial" => 0.5, @@ -151,7 +165,7 @@ function _main( data = nothing; debug = false, plot_solution = false, - Plotter = nothing, + Plotter = CairoMakie, kwargs... ) if isnothing(data) @@ -173,7 +187,6 @@ function _main( θ_spatial = data["θ_spatial"] multi_indices = Array{Array{Int, 1}, 1}(data["initial_modes"]) nsamples = data["nsamples"] - f! = data["f"] bonus_quadorder_a = data["bonus_quadorder_a"] bonus_quadorder_f = data["bonus_quadorder_f"] factor_tail = data["factor_tail"] @@ -241,9 +254,12 @@ function _main( FES = [FESpace{FETypes[1]}(xgrid), FESpace{FETypes[2]}(xgrid)] unames = ["u", "p"] + Samples, _ = sample_distribution(TensorBasis, 2; M = 2, Mweights = 2) + f! = isnothing(data["f"]) ? prepare_data(C, Samples) : data["f"] + sol = SGFEVector(FES, TensorBasis; active_modes = 1:length(multi_indices), unames = unames) - time_solve = @elapsed bdofs = solve!(problem, sol, C; rhs = (f!), bonus_quadorder_a = bonus_quadorder_a, bonus_quadorder_f = bonus_quadorder_f, use_iterative_solver = use_iterative_solver, debug = debug) + time_solve = @elapsed bdofs = solve!(problem, sol, C; (rhs!) = (f!), bonus_quadorder_a = bonus_quadorder_a, bonus_quadorder_f = bonus_quadorder_f, use_iterative_solver = use_iterative_solver) df[lvl, :time_solve] = time_solve if isnothing(Plotter) && plot_solution @@ -251,13 +267,15 @@ function _main( if debug @show extrema(view(sol[j])) end - println(stdout, unicode_scalarplot(sol[j]; title = "MI $(multi_indices[j])")) + println(stdout, unicode_scalarplot(sol[1]; title = "MI $(multi_indices[1])")) end elseif plot_solution - plot_modes(sol; Plotter = Plotter, ncols = 4) + p = plot_modes(sol; Plotter = CairoMakie, ncols = 4) + display(p) end + weightederrorH1, weightederrorL2u, weightederrorL2p, uniformerrorH1, uniformerrorL2u, uniformerrorL2p = calculate_sampling_error_2( - sol, C; problem = problem, metrics_configurations = stokes_metrics_configuration, rhs = (f!), order = order + 1, nsamples, debug + sol, C; problem, metrics_configurations = stokes_metrics_configuration, (rhs!) = (f!), order = order + 1, nsamples, debug ) tail_extension = data["tail_extension"] @@ -430,7 +448,7 @@ function produce_plots(; ylabel = L"|| u - u_h ||_{L^2}" elseif cols[1] == :exact_error_stress ylabel = L"|| u - u_h ||_A" - elseif col[1] == :exect_error_p + elseif col[1] == :exact_error_p ylabel = L"|| p - p_h ||_{L^2}" elseif cols[1] == :estimate ylabel = L"\eta" diff --git a/scripts/stokes_simple.jl b/scripts/stokes_simple.jl index 9257b2a..a04685b 100644 --- a/scripts/stokes_simple.jl +++ b/scripts/stokes_simple.jl @@ -5,64 +5,63 @@ using ExtendableFEM using ExtendableFEMBase using ExtendableGrids using GridVisualize +using Symbolics -function f!(result, qpinfo) - x = qpinfo.x[1] - y = qpinfo.x[2] - result[1] = 5(x^4) + 12(x^2) * ((-1 + x)^2) * (-1 + y) + 12(x^2) * ((-1 + x)^2) * y + 4(x^2) * y * ((-1 + y)^2) + 4(x^2) * (-1 + y) * (y^2) + 4((-1 + x)^2) * (-1 + y) * (y^2) + 4((-1 + x)^2) * y * ((-1 + y)^2) + 16x * (-1 + x) * (-1 + y) * (y^2) + 16x * (-1 + x) * y * ((-1 + y)^2) - result[2] = 5(y^4) - 4x * ((-1 + x)^2) * ((-1 + y)^2) - 4(-1 + x) * (x^2) * ((-1 + y)^2) - 4x * ((-1 + x)^2) * (y^2) - 4(-1 + x) * (x^2) * (y^2) - 16(-1 + x) * (x^2) * y * (-1 + y) - 16x * ((-1 + x)^2) * y * (-1 + y) - 12x * (y^2) * ((-1 + y)^2) - 12(-1 + x) * (y^2) * ((-1 + y)^2) - return nothing - #result[1] = 5(x^4) - 0.001(-12(x^2) * ((-1 + x)^2) * (-1 + y) - 12(x^2) * ((-1 + x)^2) * y - 4(x^2) * y * ((-1 + y)^2) - 4(x^2) * (-1 + y) * (y^2) - 4((-1 + x)^2) * (-1 + y) * (y^2) - 4((-1 + x)^2) * y * ((-1 + y)^2) - 16x * (-1 + x) * (-1 + y) * (y^2) - 16x * (-1 + x) * y * ((-1 + y)^2)) - #result[2] = 5(y^4) - 0.001(4x * ((-1 + x)^2) * ((-1 + y)^2) + 4(-1 + x) * (x^2) * ((-1 + y)^2) + 4x * ((-1 + x)^2) * (y^2) + 4(-1 + x) * (x^2) * (y^2) + 16(-1 + x) * (x^2) * y * (-1 + y) + 16x * ((-1 + x)^2) * y * (-1 + y) + 12x * (y^2) * ((-1 + y)^2) + 12(-1 + x) * (y^2) * ((-1 + y)^2)) - #result[1] = 1000 * (1 - x) * x * y * (1 - y) - #result[2] = 1000 * (1 - x) * x * y * (1 - y) +function prepare_data() + @variables x1, x2 + + u = [0, (1 / 2) * (1 - x1^2)] + eval_u! = build_function(u, x1, x2, expression = Val{false})[2] + + return (result, _) -> (result .= 0), + (result, qpinfo) -> (eval_u!(result, qpinfo.x...)) end function main(; problem = StokesProblemPrimal, - nrefs = 3, # number of uniform refinements of the initial grid - order = 2, # polynomial order of the FEspaces - decay = 2.0, # decay factor for the random coefficient - mean = 1, # mean value of coefficient - domain = "square", # domain, e.g., "square" or "lshape" - initial_modes = [[0], [1, 0], [0, 1], [2, 0], [0, 0, 1]], # initial multi-indices for stochastic basis - (f!) = (f!), # right-hand side function - use_iterative_solver = true, # use iterative solver ? (otherwise direct) - nsamples = 20, + nrefs = 3, ## number of uniform refinements of the initial grid + order = 2, ## polynomial order of the FEspaces + domain = "square", ## domain, e.g., "square" or "lshape" + modes = [[0], [1]], ## initial multi-indices for stochastic basis + use_iterative_solver = true, ## use iterative solver? (otherwise direct) + nsamples = 50, Plotter = nothing, - debug = false, ) ## prepare stochastic coefficient - C = StochasticCoefficientCosinus(; τ = 0.9, decay = decay, mean = mean) + C = StochasticCoefficientConstants() ## prepare grid - xgrid = if domain == "square" - uniform_refine(grid_unitsquare(Triangle2D), nrefs) + (xgrid, boundary_regions) = if domain == "square" + (uniform_refine(grid_unitsquare(Triangle2D), nrefs), 1:4) elseif domain == "lshape" - uniform_refine(grid_lshape(Triangle2D), nrefs) + (uniform_refine(grid_lshape(Triangle2D), nrefs), 1:8) else error("unknown domain: $domain") end ## prepare stochastic basis - multi_indices = Array{Array{Int, 1}, 1}(initial_modes) + multi_indices = Array{Array{Int, 1}, 1}(modes) prepare_multi_indices!(multi_indices) M = maximum(length.(multi_indices)) OBType = LegendrePolynomials ansatz_deg = maximum([maximum(multi_indices[k]) for k in 1:length(multi_indices)]) + 4 TensorBasis = TensorizedBasis(OBType, M, ansatz_deg, 2 * ansatz_deg, 2 * ansatz_deg, multi_indices = multi_indices) + ## prepare synthetic problem data + f!, exact_u! = prepare_data() + bonus_quadorder_f = 2 + ## prepare FE spaces FETypes = (H1BR{2}, L2P0{1}) FES = [FESpace{FETypes[1]}(xgrid), FESpace{FETypes[2]}(xgrid)] unames = ["u", "p"] ## create solution vector - sol = SGFEVector(FES, TensorBasis; active_modes = 1:length(multi_indices), unames = unames) + sol = SGFEVector(FES, TensorBasis; active_modes = 1:length(multi_indices), unames) ## solve problem @info "Solving..." - solve!(problem, sol, C; rhs = (f!), use_iterative_solver = use_iterative_solver) + solve!(problem, sol, f!, C; (exact_boundary!) = (exact_u!), bonus_quadorder_f, use_iterative_solver) ## plot solution if !isnothing(Plotter) @@ -71,12 +70,21 @@ function main(; end ## compute exact error (by MC sampling) - weightederrorH1, weightederrorL2u, weightederrorL2p, uniformerrorH1, uniformerrorL2u, uniformerrorL2p = calculate_sampling_error_2( - sol, C; problem = problem, metrics_configurations = stokes_metrics_configuration, rhs = (f!), order = order + 1, nsamples, debug + weightederrorH1, weightederrorL2u, weightederrorL2p, uniformerrorH1, uniformerrorL2u, uniformerrorL2p = + calculate_sampling_error_2( + sol, + f!, + C; + problem, + metrics_configurations = stokes_metrics_configuration, + (exact_boundary!) = (exact_u!), + boundary_regions, + order = order + 1, + nsamples, ) @info "RESULTS - || ∇(u-u_h) || (w,u) = $(sqrt(weightederrorH1[end])), $(sqrt(uniformerrorH1[end])) + || ∇(u - u_h) || (w,u) = $(sqrt(weightederrorH1[end])), $(sqrt(uniformerrorH1[end])) || u - u_h || (w,u) = $(sqrt(weightederrorL2u[end])), $(sqrt(uniformerrorL2u[end])) || p - p_h || (w,u) = $(sqrt(weightederrorL2p[end])), $(sqrt(uniformerrorL2p[end]))" diff --git a/src/ExtendableASGFEM.jl b/src/ExtendableASGFEM.jl index 5719f2d..ab002c8 100644 --- a/src/ExtendableASGFEM.jl +++ b/src/ExtendableASGFEM.jl @@ -7,12 +7,12 @@ using ExtendableFEM: ExtendableFEM, BilinearOperator, FaceInterpolator, HomogeneousBoundaryData, ItemIntegrator, ItemIntegratorDG, L2NormIntegrator, LinearOperator, ProblemDescription, Reconstruct, Unknown, ZeroMeanValueRestriction, apply, assemble!, assign_operator!, assign_restriction!, - assign_unknown!, div, grad, id, jump, plot + assign_unknown!, div, grad, id, jump, plot, InterpolateBoundaryData using ExtendableFEMBase: ExtendableFEMBase, BFaceDofs, CellDofs, Divergence, FEEvaluator, FEMatrix, FESpace, FEVector, FEVectorBlock, Gradient, H1P1, H1Pk, H1P2B, L2P1, HDIVRT0, HDIVRT1, HDIVRTk, H1BR, Identity, L2P0, Laplacian, QuadratureRule, _addnz, addblock!, - addblock_matmul!, eval_febe!, fill!, get_ndofs, + addblock_matmul!, eval_febe!, fill!, get_ndofs, ON_BFACES, get_polynomialorder, unicode_scalarplot, update_basis! using ExtendableGrids: ExtendableGrids, Adjacency, BFaceFaces, CellFaces, @@ -25,7 +25,8 @@ using ExtendableGrids: ExtendableGrids, Adjacency, BFaceFaces, CellFaces, max_num_targets_per_source, num_cells, num_nodes, num_sources, num_targets, unique, update_trafo! using ExtendableSparse: ExtendableSparse, ExtendableSparseMatrix, flush! -using GridVisualize: GridVisualize, GridVisualizer, scalarplot, scalarplot! +using GridVisualize: GridVisualize, GridVisualizer, scalarplot, scalarplot!, + GridVisualize using IterativeSolvers: IterativeSolvers using Krylov: Krylov using LinearAlgebra: LinearAlgebra, SymTridiagonal, dot, eigvals, eigvecs, @@ -70,6 +71,7 @@ export set_sample! include("coefficients/coefficients.jl") export AbstractStochasticCoefficient export StochasticCoefficientCosinus, SingleStochasticCoefficient +export StochasticCoefficientConstants export get_am!, get_gradam! export get_a!, get_expa! export get_am_x, get_gradam_x_sigma, get_gradam_x_u diff --git a/src/coefficients/coefficients.jl b/src/coefficients/coefficients.jl index 908c35d..4eb9c52 100644 --- a/src/coefficients/coefficients.jl +++ b/src/coefficients/coefficients.jl @@ -14,7 +14,7 @@ of ``a``). """ abstract type AbstractStochasticCoefficient{T} end include("cosinus.jl") - +include("constant_coefficients.jl") """ $(TYPEDSIGNATURES) diff --git a/src/coefficients/constant_coefficients.jl b/src/coefficients/constant_coefficients.jl new file mode 100644 index 0000000..856df84 --- /dev/null +++ b/src/coefficients/constant_coefficients.jl @@ -0,0 +1,46 @@ +""" +$(TYPEDEF) + +Expansion of the form + + a(x,y) = a_0(x) + ∑_m y_m a_m(x) + +where the `a_m` are of cosinus type. + +""" +struct StochasticCoefficientConstants{T} <: AbstractStochasticCoefficient{T} + constants::Vector{T} +end + +maxm(SC::StochasticCoefficientConstants) = length(SC.constants) +meanvalue(SC::StochasticCoefficientConstants) = SC.constants[1] + + +""" +$(TYPEDSIGNATURES) + +constructor for StochasticCoefficientConstants of type `T` (default = Float64), +where `decay` (default = 2) steers the decay of the coefficient basis functions (the larger the faster), +`mean` (default = 0) is the mean value of the coefficient, `maxm` (default = 100) is the maximal number of stochastic random variables, +and `τ` (default = 1) is a uniform scaling factors + +""" +function StochasticCoefficientConstants(; constants = [1.0, 0.2]) + return StochasticCoefficientConstants{eltype(constants)}(constants) +end + + +function get_am!(result, x, m, SC::StochasticCoefficientConstants) + result[1] = m < length(SC.constants) ? SC.constants[m + 1] : 0 + return nothing +end + +function get_gradam!(result, x, m, SC::StochasticCoefficientConstants) + fill!(result, 0) + return nothing +end + +function Base.show(io::IO, SC::StochasticCoefficientConstants) + println(io, "COEFFICIENT DATA") + return println(io, "constants = $(SC.constants)") +end diff --git a/src/estimate.jl b/src/estimate.jl index 48197d4..fec7b0e 100644 --- a/src/estimate.jl +++ b/src/estimate.jl @@ -381,6 +381,7 @@ end # TODO: Implement StokesProblemPrimal a posteriori estimator # This is a dummy estimator which enforces repeated spatial refinement +level = 0 function estimate(::Type{StokesProblemPrimal}, sol::SGFEVector, C::AbstractStochasticCoefficient; rhs = nothing, bonus_quadorder = 1, tail_extension = 5) FES = sol.FES_space[1] xgrid = FES.xgrid @@ -401,9 +402,10 @@ function estimate(::Type{StokesProblemPrimal}, sol::SGFEVector, C::AbstractStoch actives = union(active_int, active_bnd) for j in 1:length(nmodes_extended) if j in actives - eta4modes[j] = 10000 + eta4modes[j] = ((level % 2) == 0) ? 10000 : 0 end end + global level += 1 return eta4modes, eta4cell, multi_indices_extended end diff --git a/src/modelproblems/solvers_stokes_primal.jl b/src/modelproblems/solvers_stokes_primal.jl index 1fd648f..dd1220c 100644 --- a/src/modelproblems/solvers_stokes_primal.jl +++ b/src/modelproblems/solvers_stokes_primal.jl @@ -22,11 +22,6 @@ struct StokesPrimalPreconditioner{Tv, FAC} end function stokesPrimalPreconditioner(A0::ExtendableSparseMatrix{Tv, Ti}, B::ExtendableSparseMatrix{Tv, Ti}, bdofs, nmodes, vsize) where {Tv, Ti} - #for dof in bdofs - # A0[dof, dof] = 1.0e60 - #end - #flush!(A0) - DA::Array{Tv, 1} = zeros(Tv, size(A0, 1)) for j in 1:length(DA) DA[j] = A0[j, j] @@ -52,7 +47,6 @@ function stokesPrimalPreconditioner(A0::ExtendableSparseMatrix{Tv, Ti}, B::Exten # compute LU factorisation of S flush!(S) LUS = lu(S.cscmatrix) - #LUA = lu(A0.cscmatrix) # temporary storage array for solver temp = zeros(Tv, size(B, 2)) @@ -67,6 +61,7 @@ end DA::Array{Tv, 1} = C.DA temp::Array{Tv, 1} = C.temp nmodes::Int = C.nmodes + bdofs = C.bdofs vsize::Array{Int, 1} = C.vsize for mu in 1:nmodes # upper left block of preconditioner (I ⊗ A_diag) @@ -83,6 +78,11 @@ end ldiv!(temp, C.LUS, view(b, a:c)) y[a:c] .= temp end + + a = (mu - 1) * vsize[1] + 1 + for dof in bdofs + y[a + dof - 1] = 0 + end end return y @@ -138,8 +138,9 @@ function LinearAlgebra.mul!(Ax::Vector{Tv}, S::StokesPrimal{Tv, MT, VT, GT}, x) end for dof in bdofs - Ax[a + dof - 1] = 0 + Ax[a + dof - 1] = (mu == 1) ? 1.0e60 : 0 end + Ax[nmodes * vsize[1] + 1] = 0 # Wegen pressure end return nothing @@ -148,18 +149,13 @@ end Base.eltype(S::StokesPrimal) = typeof(S).parameters[1] Base.size(S::StokesPrimal) = S.nmodes .* (size(S.A0.entries) .+ size(S.B.entries)[2]) -function solve_stokes_primal!(SolutionSGFEM::SGFEVector, A0, A, B, b0, G, nmodes, bfac; atol = 1.0e-14, rtol = 1.0e-14) +function solve_stokes_primal!( + SolutionSGFEM::SGFEVector, A0, A, B, b0, G, nmodes, bdofs; + atol = 1.0e-14, rtol = 1.0e-14 + ) ## create fullmatrix-free matrix evaluator @info "Solving StochasticFEM iteratively and matrix-free (ndofs = $(length(SolutionSGFEM.entries)))..." vsize = [SolutionSGFEM[1].FES.ndofs, SolutionSGFEM[nmodes + 1].FES.ndofs] - ## boundary data - bfacedofs = SolutionSGFEM.FES_space[1][BFaceDofs] - nbfaces = num_sources(bfacedofs) - bdofs = [] - for bface in 1:nbfaces - append!(bdofs, view(bfacedofs, :, bface)) - end - bdofs = unique(bdofs) S = StokesPrimal{eltype(G), typeof(A0), typeof(b0), typeof(G)}(A0, A, B, G, bdofs, nmodes, vsize) @info "...initializing Preconditioner" @@ -167,42 +163,46 @@ function solve_stokes_primal!(SolutionSGFEM::SGFEVector, A0, A, B, b0, G, nmodes ## right-hand side b = deepcopy(SolutionSGFEM) - fill!(b.entries, 0) - addblock!(b[1], b0[1]; factor = bfac) - for m in 1:nmodes - for dof in bdofs - b[m][dof] = 0 - end + addblock!(b[1], b0[1]) + for dof in bdofs + b[1][dof] *= 1.0e60 end ## solve + addblock!(SolutionSGFEM[1], b0[1]) @info "...starting right-conditioned GMRES" x, history = Krylov.gmres(S, b.entries, SolutionSGFEM.entries; ldiv = true, atol = atol, rtol = rtol, M = P) SolutionSGFEM.entries .= x @show history - ## Pressure mean von pressure modes abziehen, damit Vorfaktor wie bei use_iterative_solver=false ist. - ## https://wias-pdelib.github.io/ExtendableFEM.jl/stable/module_examples/Example252_NSEPlanarLatticeFlow/ + # enforce uniqueness of the pressure solution xgrid = SolutionSGFEM[1].FES.xgrid for i in 1:nmodes pintegrate = ItemIntegrator([id(1)]) - pmean = sum(ExtendableFEM.evaluate(pintegrate, [SolutionSGFEM.FEVectorBlocks[nmodes + i]])) / sum(xgrid[CellVolumes]) + pmean = sum(ExtendableFEM.evaluate(pintegrate, [SolutionSGFEM.FEVectorBlocks[nmodes + i]])) / + sum(xgrid[CellVolumes]) view(SolutionSGFEM.FEVectorBlocks[nmodes + i]) .-= pmean end ## check residual + FES = SolutionSGFEM.FES_space Ax = zero(SolutionSGFEM.entries) mul!(Ax, S, SolutionSGFEM.entries) - @info "solver residual = $(sqrt(sum((Ax - b.entries) .^ 2)))" + residual = Ax - b.entries + a = FES[1].ndofs + b = FES[2].ndofs + for m in 1:nmodes + residual[a * (m - 1) .+ bdofs] .= 0 + residual[a * nmodes + b * (m - 1) + 1] = 0 + end + @info "linear residual = $(sqrt(sum(residual .^ 2)))" - return bdofs + return nothing end -function solve_stokes_primal_full!(SolutionSGFEM::SGFEVector, A0, A, B, b0, G, nmodes, rhsfac) - M::Int = length(A) # size(G,1) / nmodes +function solve_stokes_primal_full!(SolutionSGFEM::SGFEVector, A0, A, B, b0, G, nmodes, bdofs) + M::Int = length(A) FES = SolutionSGFEM.FES_space - multi_indices = SolutionSGFEM.TB.multi_indices - nmodes = num_multiindices(SolutionSGFEM) bigFES = Array{FESpace{Float64, Int32}, 1}([FES[1] for j in 1:nmodes]) append!(bigFES, [FES[2] for j in 1:nmodes]) @@ -229,38 +229,43 @@ function solve_stokes_primal_full!(SolutionSGFEM::SGFEVector, A0, A, B, b0, G, n end flush!(bigS.entries) - ## boundary data - bfacedofs = FES[1][BFaceDofs] - nbfaces = num_sources(bfacedofs) - bdofs = [] - for bface in 1:nbfaces - append!(bdofs, view(bfacedofs, :, bface)) + ## boundary condition for bigS + for dof in bdofs + bigb[1][dof] *= 1.0e60 + for m in 1:nmodes + bigS[m, m][dof, dof] = 1.0e60 + end end - unique!(bdofs) + flush!(bigS.entries) + ## include integral mean condition for each pressure block for m in 1:nmodes - for dof in bdofs - bigS[m, m][dof, dof] = 1.0e60 - bigb[m][dof] = 0 - end + bigS[nmodes + m, nmodes + m][1, 1] = 1.0e60 + bigb[nmodes + m][1] = 0 end flush!(bigS.entries) @info "Solving StochasticFEM with full matrix..." SolutionSGFEM.entries .= bigS.entries \ bigb.entries + ## enforce uniqueness of the pressure solution xgrid = SolutionSGFEM[1].FES.xgrid for i in 1:nmodes - pintegrate = ItemIntegrator([id(1)]) - pmean = sum(ExtendableFEM.evaluate(pintegrate, [SolutionSGFEM.FEVectorBlocks[nmodes + i]])) / sum(xgrid[CellVolumes]) + p_integrator = ItemIntegrator([id(1)]) + pmean = sum(ExtendableFEM.evaluate(p_integrator, [SolutionSGFEM.FEVectorBlocks[nmodes + i]])) / + sum(xgrid[CellVolumes]) view(SolutionSGFEM.FEVectorBlocks[nmodes + i]) .-= pmean end + ## calculate residual residual = bigS.entries * SolutionSGFEM.entries .- bigb.entries + a = FES[1].ndofs + b = FES[2].ndofs for m in 1:nmodes - residual[FES[1].ndofs * (m - 1) .+ bdofs] .= 0 + residual[a * (m - 1) .+ bdofs] .= 0 + residual[a * nmodes + b * (m - 1) + 1] = 0 end - println("linear residual = $(sqrt(sum(residual .^ 2)))") + @info "linear residual = $(sqrt(sum(residual .^ 2)))" - return bdofs + return nothing end diff --git a/src/modelproblems/stokes_primal.jl b/src/modelproblems/stokes_primal.jl index 1e635ac..74c7487 100644 --- a/src/modelproblems/stokes_primal.jl +++ b/src/modelproblems/stokes_primal.jl @@ -11,11 +11,21 @@ abstract type StokesProblemPrimal <: AbstractModelProblem end include("solvers_stokes_primal.jl") ## deterministic problem description -function deterministic_problem(::Type{StokesProblemPrimal}, C::AbstractStochasticCoefficient, sample_pointer; (get_a!) = (get_a!), rhs = nothing, bonus_quadorder_a = 2, bonus_quadorder_f = 0) +function deterministic_problem( + ::Type{StokesProblemPrimal}, + rhs!::Function, + C::AbstractStochasticCoefficient, + sample_pointer; + (get_a!) = (get_a!), + (exact_boundary!) = nothing, + reconstruct = true, + boundary_regions = 1:4, + bonus_quadorder_a = 2, + bonus_quadorder_f = 0, + ) get_ν! = get_a!(C) - + ν = zeros(1) function stokes_kernel!(result, input, qpinfo) - ν = zeros(1) get_ν!(ν, qpinfo.x, sample_pointer) result[1] = ν[1] * input[1] - input[5] result[2] = ν[1] * input[2] @@ -31,11 +41,25 @@ function deterministic_problem(::Type{StokesProblemPrimal}, C::AbstractStochasti assign_unknown!(PD, u) assign_unknown!(PD, p) assign_operator!(PD, BilinearOperator(stokes_kernel!, [grad(u), id(p)]; bonus_quadorder = bonus_quadorder_a)) - if rhs !== nothing - # TODO: Define reconstruct::Bool switch for reconstruction operator - assign_operator!(PD, LinearOperator(rhs, [apply(u, Reconstruct{HDIVRT1{2}, Identity})]; bonus_quadorder = bonus_quadorder_f)) + assign_operator!( + PD, + LinearOperator( + rhs!, [reconstruct ? apply(u, Reconstruct{HDIVRT1{2}, Identity}) : id(u)]; + bonus_quadorder = bonus_quadorder_f + ) + ) + # Is a non-zero boundary condition or some non-standard definition required? + if !isnothing(exact_boundary!) + assign_operator!( + PD, + InterpolateBoundaryData( + u, exact_boundary!; + regions = boundary_regions, bonus_quadorder = bonus_quadorder_f + ) + ) + else # Define zero boundary if no additional constraints are defined + assign_operator!(PD, HomogeneousBoundaryData(u; regions = boundary_regions)) end - assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:4)) assign_restriction!(PD, ZeroMeanValueRestriction(p)) return PD, u @@ -44,13 +68,14 @@ end ## solver for stochastic Galerkin problem function solve!( ::Type{StokesProblemPrimal}, - sol::SGFEVector, ## target SGFEM vector - C::AbstractStochasticCoefficient; ## stochastic coefficient a - bonus_quadorder_a = 2, ## additional quadrature order for grad(am) - bonus_quadorder_f = 0, ## additional quadrature order for rhs f - rhs = nothing, - debug = false, - use_iterative_solver = true + sol::SGFEVector, ## target SGFEM vector + rhs!::Function, ## kernel for right-hand side calculation + C::AbstractStochasticCoefficient; ## stochastic coefficient ν + (exact_boundary!) = nothing, ## kernel function for boundary values + bonus_quadorder_a = 2, ## additional quadrature order for grad(am) + bonus_quadorder_f = 0, ## additional quadrature order for rhs f + use_iterative_solver = true, ## use iterative solver? (otherwise direct) + reconstruct = false, ## apply HDIVRT0{2} reconstruction operator to rhs? ) FES = sol.FES_space TB = sol.TB @@ -68,14 +93,30 @@ function solve!( ## Deterministic B = (∇ ⋅ v, q) term B = FEMatrix(FES[1], FES[2]) - assemble!(B, BilinearOperator([div(1)], [id(1)]; factor = -1)) # Checken! + assemble!(B, BilinearOperator([div(1)], [id(1)]; factor = -1)) ## assemble right-hand side (f, v) - if rhs !== nothing - b0 = FEVector(FES[1]) - assemble!(b0, LinearOperator(rhs, [apply(1, Reconstruct{HDIVRT0{2}, Identity})]; bonus_quadorder = bonus_quadorder_f)) - else - @error "need right-hand side" + b0 = FEVector(FES[1]) + assemble!( + b0, + LinearOperator( + rhs!, [reconstruct ? apply(1, Reconstruct{HDIVRT0{2}, Identity}) : id(1)]; + bonus_quadorder = bonus_quadorder_f + ) + ) + + ## calculate bdofs for boundary conditions and future use + bfacedofs = FES[1][BFaceDofs] + nbfaces = num_sources(bfacedofs) + bdofs = [] + for bface in 1:nbfaces + append!(bdofs, view(bfacedofs, :, bface)) + end + unique!(bdofs) + + ## set appropriate boundary on the right-hand side + if !isnothing(exact_boundary!) + interpolate!(b0[1], ON_BFACES, exact_boundary!; bonus_quadorder = bonus_quadorder_f) end ## get G matrix @@ -83,11 +124,11 @@ function solve!( nmodes = num_multiindices(TB) ## solve - if use_iterative_solver - @time bdofs = solve_stokes_primal!(sol, A0, A, B, b0, G, nmodes, 1) - else - @time bdofs = solve_stokes_primal_full!(sol, A0, A, B, b0, G, nmodes, 1) - end + @time ( + use_iterative_solver + ? solve_stokes_primal! + : solve_stokes_primal_full! + )(sol, A0, A, B, b0, G, nmodes, bdofs) return bdofs end diff --git a/src/plots.jl b/src/plots.jl index ee6a636..86a0fac 100644 --- a/src/plots.jl +++ b/src/plots.jl @@ -58,5 +58,6 @@ function plot_basis(ONB::ONBasis{T, OBT, npoly, nquad}; Plotter = nothing, resol for j in 1:npoly scalarplot!(p[1, 1], xgrid, vals4xref; clear = false, Plotter = Plotter, label = "$j", xlimits = (-1, 1), limits = (-4, 4)) end + GridVisualize.reveal(p) return end diff --git a/src/sampling_error.jl b/src/sampling_error.jl index 444d6b4..1811e3e 100644 --- a/src/sampling_error.jl +++ b/src/sampling_error.jl @@ -1,92 +1,45 @@ -## kernel for L2 error of stress, i.e || grad(u) - grad(u_h) || -function data_error_stress(::Type{<:PoissonProblemPrimal}, dim, C::AbstractStochasticCoefficient, sample_pointer) - function closure(result, input, qpinfo) - result[1] = 1 - if dim == 1 - result[2] = (input[1] - input[2])^2 - elseif dim == 2 - result[2] = (input[1] - input[3])^2 + (input[2] - input[4])^2 - end - result[1] *= result[2] - return nothing +using PythonPlot + +function create_squared_l2_dist_closure!(dim = nothing) + N = !isnothing(dim) ? dim : trunc(Int, length(input) / 2) + return function (result, input, _) + return result[1] = sum((input[1:N] - input[(N + 1):(N << 1)]) .^ 2) end - return closure, [grad(1), grad(2)], [(1, 1), (2, 1)] end -function data_error_stress(::Type{<:LogTransformedPoissonProblemDual}, dim, C::AbstractStochasticCoefficient, sample_pointer) - function closure(result, input, qpinfo) - result[1] = 1 - if dim == 1 - result[2] = (input[1] - input[2])^2 - elseif dim == 2 - result[2] = (input[1] - input[3])^2 + (input[2] - input[4])^2 - end - result[1] *= result[2] - return nothing - end - return closure, [id(1), id(2)], [(1, 1), (2, 1)] +## kernel for L2 error of stress, i.e || grad(u) - grad(u_h) || +function data_error_stress(::Union{Type{<:PoissonProblemPrimal}, Type{<:LogTransformedPoissonProblemPrimal}, Type{<:StokesProblemPrimal}}, dim, C::AbstractStochasticCoefficient, sample_pointer) + return create_squared_l2_dist_closure!(dim), [grad(1), grad(2)], [(1, 1), (2, 1)] end -function data_error_stress(::Union{Type{<:StokesProblemPrimal}, Type{<:LogTransformedPoissonProblemPrimal}}, dim, C::AbstractStochasticCoefficient, sample_pointer) - function closure(result, input, qpinfo) - result[1] = 1 - if dim == 1 - result[2] = (input[1] - input[2])^2 - elseif dim == 2 - result[2] = (input[1] - input[3])^2 + (input[2] - input[4])^2 - end - result[1] *= result[2] - return nothing - end - return closure, [grad(1), grad(2)], [(1, 1), (2, 1)] +function data_error_stress(::Type{<:LogTransformedPoissonProblemDual}, dim, C::AbstractStochasticCoefficient, sample_pointer) + return create_squared_l2_dist_closure!(dim), [id(1), id(2)], [(1, 1), (2, 1)] end ## kernel for L2 error || u - u_h || -function data_error_u(::Union{Type{<:PoissonProblemPrimal}, Type{<:LogTransformedPoissonProblemPrimal}}, dim, C::AbstractStochasticCoefficient, sample_pointer) - function closure(result, input, qpinfo) - result[1] = (input[1] - input[2])^2 - return nothing - end - return closure, [id(1), id(2)], [(1, 1), (2, 1)] +function data_error_u(::Union{Type{<:PoissonProblemPrimal}, Type{<:LogTransformedPoissonProblemPrimal}, Type{<:StokesProblemPrimal}}, dim, C::AbstractStochasticCoefficient, sample_pointer) + return create_squared_l2_dist_closure!(dim), [id(1), id(2)], [(1, 1), (2, 1)] end function data_error_u(::Type{<:LogTransformedPoissonProblemDual}, dim, C::AbstractStochasticCoefficient, sample_pointer) - function closure(result, input, qpinfo) - result[1] = (input[1] - input[2])^2 - return nothing - end - return closure, [id(1), id(2)], [(1, 2), (2, 1)] -end - -function data_error_u(::Type{<:StokesProblemPrimal}, dim, C::AbstractStochasticCoefficient, sample_pointer) - function closure(result, input, qpinfo) - result[1] = (input[1] - input[2])^2 - return nothing - end - return closure, [id(1), id(2)], [(1, 1), (2, 1)] + return create_squared_l2_dist_closure!(dim), [id(1), id(2)], [(1, 2), (2, 1)] end +## kernel for L2 error || p - p_h || function data_error_p(::Type{<:StokesProblemPrimal}, dim, C::AbstractStochasticCoefficient, sample_pointer) - function closure(result, input, qpinfo) - result[1] = (input[1] - input[2])^2 - return nothing - end - return closure, [id(1), id(2)], [(1, 2), (2, 2)] + return create_squared_l2_dist_closure!(dim), [id(1), id(2)], [(1, 2), (2, 2)] end -FES4sampling(::Type{PoissonProblemPrimal}, dim, xgrid, order) = [FESpace{H1Pk{1, dim, order}}(xgrid)] -FES4sampling(::Type{LogTransformedPoissonProblemPrimal}, dim, xgrid, order) = [FESpace{H1Pk{1, dim, order}}(xgrid)] -FES4sampling(::Type{LogTransformedPoissonProblemDual}, dim, xgrid, order) = [FESpace{HDIVRTk{dim, order}}(xgrid), FESpace{H1Pk{1, dim, order}}(xgrid; broken = true)] -FES4sampling(::Type{StokesProblemPrimal}, _, xgrid, order) = [FESpace{H1P2B{2, 2}}(xgrid), FESpace{L2P1{1}}(xgrid)] - u_with_stress_metric_configuration = [ Dict( "name" => "L2stress", "closure" => data_error_stress, + "dim" => 4, ), Dict( "name" => "L2u", "closure" => data_error_u, + "dim" => 2, ), ] @@ -94,17 +47,26 @@ stokes_metrics_configuration = [ Dict( "name" => "L2stress", "closure" => data_error_stress, + "dim" => 4, ), Dict( "name" => "L2u", "closure" => data_error_u, + "dim" => 2, ), Dict( "name" => "L2p", "closure" => data_error_p, + "dim" => 1, ), ] + +FES4sampling(::Type{PoissonProblemPrimal}, dim, xgrid, order) = [FESpace{H1Pk{1, dim, order}}(xgrid)] +FES4sampling(::Type{LogTransformedPoissonProblemPrimal}, dim, xgrid, order) = [FESpace{H1Pk{1, dim, order}}(xgrid)] +FES4sampling(::Type{LogTransformedPoissonProblemDual}, dim, xgrid, order) = [FESpace{HDIVRTk{dim, order}}(xgrid), FESpace{H1Pk{1, dim, order}}(xgrid; broken = true)] +FES4sampling(::Type{StokesProblemPrimal}, _, xgrid, _) = [FESpace{H1P2B{2, 2}}(xgrid), FESpace{L2P1{1}}(xgrid)] + """ $(TYPEDSIGNATURES) @@ -142,7 +104,7 @@ function calculate_sampling_error( SolutionSGFEM::SGFEVector, C::AbstractStochasticCoefficient; problem = LogTransformedPoissonProblemPrimal, - bonus_quadorder_a = 10, + bonus_quadorder_a = 2, bonus_quadorder_f = 0, order = 2, rhs = nothing, @@ -277,14 +239,16 @@ end function calculate_sampling_error_2( SolutionSGFEM::SGFEVector, + rhs!::Function, C::AbstractStochasticCoefficient; metrics_configurations = u_with_stress_metric_configuration, problem = LogTransformedPoissonProblemPrimal, - bonus_quadorder_a = 10, + (exact_boundary!) = nothing, + boundary_regions = 1:4, + dim = size(SolutionSGFEM.FES_space[1].xgrid[Coordinates], 1), + bonus_quadorder_a = 2, bonus_quadorder_f = 0, order = 2, - rhs = nothing, - dim = size(SolutionSGFEM.FES_space[1].xgrid[Coordinates], 1), Msamples = maxm(C), parallel_sampling = true, dimensionwise_error = true, @@ -299,9 +263,7 @@ function calculate_sampling_error_2( xgrid = FES[1].xgrid sol_sgfem = SolutionSGFEM.FEV TensorBasis = SolutionSGFEM.TB - multi_indices::Array{Array{Int, 1}, 1} = TensorBasis.multi_indices M::Int = maxlength_multiindices(TensorBasis) - nmodes = num_multiindices(TensorBasis) ## generate samples Samples, weights = sample_distribution(TensorBasis, nsamples; M = Msamples, Mweights = Msamples) @@ -311,9 +273,17 @@ function calculate_sampling_error_2( ## compute deterministic solutions (in parallel) Threads.@threads for s in 1:nsamples - ## deterministic problem description - PD, u = deterministic_problem(problem, C, Samples[:, s]; rhs = rhs, bonus_quadorder_a = bonus_quadorder_a, bonus_quadorder_f = bonus_quadorder_f) + PD, _ = deterministic_problem( + problem, + rhs!, + C, + Samples[:, s]; + exact_boundary!, + boundary_regions, + bonus_quadorder_a, + bonus_quadorder_f + ) ## solve problem for the current sample FESSampling = FES4sampling(problem, dim, xgrid, order) @@ -322,6 +292,8 @@ function calculate_sampling_error_2( end println(" MC samples solved") + ExtendableFEM.plot([id(1), id(2)], sol_det[end]; Plotter = PythonPlot) + ## compute errors for each sample (sequentially) M0 = dimensionwise_error ? 0 : M error_per_run = Dict() @@ -329,13 +301,16 @@ function calculate_sampling_error_2( csample = Samples[:, 1] metric_names = [configuration["name"] for configuration in metrics_configurations] for configuration in metrics_configurations - kernel!, input, ids = configuration["closure"](problem, dim, C, csample) + kernel!, input, ids = configuration["closure"](problem, configuration["dim"], C, csample) push!( metrics, Dict( "name" => configuration["name"], "ids" => ids, - "integrator" => ItemIntegrator(kernel!, input; resultdim = dim, quadorder = 2 * order), + "integrator" => ItemIntegrator( + kernel!, input; + resultdim = configuration["dim"], quadorder = 2 * order + 1 + ), "solution" => Array{FEVectorBlock, 1}(undef, 2), ), ) @@ -376,8 +351,9 @@ function calculate_sampling_error_2( error[name_weighted] = zeros(Float64, M + 1) error[name_uniform] = zeros(Float64, M + 1) for s in 1:nsamples - error[name_uniform] .+= view(error_per_run[name], :, s) - error[name_weighted] .+= view(error_per_run[name], :, s) * weights[s] + slice = view(error_per_run[name], :, s) + error[name_uniform] .+= slice + error[name_weighted] .+= slice * weights[s] end for m in 1:(M + 1) error[name_uniform][m] /= nsamples @@ -385,5 +361,6 @@ function calculate_sampling_error_2( end end - return [error["$(name)_weighted"] for name in metric_names]..., [error["$(name)_uniform"] for name in metric_names]... + return [error["$(name)_weighted"] for name in metric_names]..., + [error["$(name)_uniform"] for name in metric_names]... end