improved documentation

Author: chmerdon

Project.toml

@@ -1,56 +1,66 @@
 name = "ExtendableASGFEM"
 uuid = "a4750d0b-645d-4217-869c-e25952806e24"
+version = "1.0.0"
 authors = ["Christian Merdon <merdon@wias-berlin.de>", "Martin Eigel <martin.eigel@wias-berlin.de>"]
-version = "1.0"
 
 [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"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DrWatson = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1"
 ExtendableFEM = "a722555e-65e0-4074-a036-ca7ce79a4aed"
 ExtendableFEMBase = "12fb9182-3d4c-4424-8fd1-727a0899810c"
 ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
 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"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 
 [compat]
-ExplicitImports = "1"
 Aqua = "0.8"
+CairoMakie = "0.13.10"
+DataFrames = "1.8.1"
 Distributed = "1.9"
 Distributions = "0.25.120"
 DocStringExtensions = "0.9.4"
+Documenter = "1.15.0"
 DoubleFloats = "1.4.3"
-Krylov = "0.10.1"
-LinearAlgebra = "1.9"
+DrWatson = "2.19.1"
+ExplicitImports = "1"
 ExtendableFEM = "1.3"
 ExtendableFEMBase = "1.3"
 ExtendableGrids = "1.13"
 ExtendableSparse = "1.7"
 GridVisualize = "1.12"
 IterativeSolvers = "0.9.2"
+Krylov = "0.10.1"
+LaTeXStrings = "1.4.0"
+LinearAlgebra = "1.9"
 Printf = "1.9"
 Random = "1.9"
 SparseArrays = "1.9"
 SpecialFunctions = "2.5.1"
 StaticArrays = "1.9.13"
 Test = "1.9"
+UnicodePlots = "3.8.1"
 julia = "^1.9"
 
-
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
-
 
 [targets]
 test = ["Aqua", "ExplicitImports", "DoubleFloats", "Test"]

docs/make.jl

@@ -28,6 +28,9 @@ makedocs(
         "Plotting" => [
             "plots.md",
         ]
+        "Scripts" => [
+            "poisson_script.md",
+        ]
     ],
 )
 

docs/src/index.md

@@ -4,6 +4,18 @@ This package provides an implementation of the stochastic Galerkin finite elemen
 
 Spatial discretization is performed using the finite element packages [ExtendableFEM.jl](https://github.com/WIAS-PDELib/ExtendableFEM.jl) and [ExtendableFEMBase.jl](https://github.com/WIAS-PDELib/ExtendableFEMBase.jl).
 
+
+## Example Script
+
+A script for running adaptive SGFEM experiments is provided in `scripts/poisson.jl`. It supports:
+- Poisson problem with linear coefficient
+- Log-transformed Poisson problem with exponential coefficient
+- Dual formulation of the log-transformed problem (WIP)
+
+The script handles both spatial and stochastic adaptivity, and includes tools for result visualization and error analysis. For a detailed description of all available parameters and usage instructions, see [Script Documentation](poisson_script.md).
+
+
+
 ## References
 
 - [1]   "Adaptive stochastic Galerkin FEM"  

docs/src/poisson_script.md

@@ -0,0 +1,69 @@
+# Script on Poisson Problems
+
+The script `scripts/poisson.jl` provides an implementation of adaptive stochastic Galerkin FEM for various formulations of the parametric Poisson problem.
+
+## Available Problem Types
+
+- `PoissonProblemPrimal`: Standard Poisson problem with linear coefficient `a`
+- `LogTransformedPoissonProblemPrimal`: Log-transformed Poisson problem with exponential coefficient `exp(a)`
+- `LogTransformedPoissonProblemDual`: Dual formulation of the log-transformed problem (WIP)
+
+## Usage
+
+The main function to run experiments is:
+```julia
+run(; problem = problem_type, kwargs...)
+```
+
+Results can be loaded and visualized using:
+```julia
+show_results(; kwargs...)    # show results and statistics
+produce_plots(; kwargs...)   # generate convergence and other plots
+```
+
+## Parameters
+
+The following parameters can be specified (with their defaults):
+
+### Problem Configuration
+- `problem = LogTransformedPoissonProblemPrimal`: Problem type to solve
+- `domain = "square"`: Domain shape ("square" or "lshape")
+- `order = 1`: Polynomial order of finite elements
+- `initial_refs = 1`: Initial uniform mesh refinements
+
+### Stochastic Parameters
+- `C = StochasticCoefficientCosinus`: Type of stochastic coefficient
+- `decay = 2`: Decay rate of the coefficient expansion
+- `mean = 0`: Mean value of the coefficient
+- `maxm = 150`: Maximum number of terms in coefficient expansion
+- `initial_modes = [[0]]`: Initial active stochastic modes
+
+### Adaptive Parameters
+- `θ_stochastic = 0.5`: Marking parameter for stochastic refinement
+- `θ_spatial = 0.5`: Marking parameter for spatial refinement
+- `factor_tail = 1`: Factor for tail estimator comparison
+- `tail_extension = [10, 2]`: Extension parameters for [0] mode and others
+- `maxdofs = 1.0e4`: Maximum degrees of freedom
+- `nsamples = 150`: Number of samples for error computation
+- `use_equilibration_estimator = false`: Use equilibration error estimator instead of standard residual estimator
+
+### Solver Options
+- `use_iterative_solver = true`: Use iterative solver (instead of a direct solver with full matrix)
+- `bonus_quadorder_a = 2`: Additional quadrature order for coefficient
+- `bonus_quadorder_f = 0`: Additional quadrature order for right-hand side
+
+### Right-Hand Side
+- `f = (result, qpinfo) -> (result[1] = 1)`: Right-hand side function
+
+### Visualization
+- `Plotter = nothing`: Plotting backend (e.g., CairoMakie)
+
+## Results Storage
+
+Results are automatically stored using [DrWatson.jl](https://github.com/JuliaDynamics/DrWatson.jl) in the `data` folder. The filename includes most problem parameters for reproducibility:
+
+```
+data/[problem]/[domain]/order=[order]_maxdofs=[maxdofs]_decay=[decay]_mean=[mean]_θ=([θ_spatial],[θ_stochastic])_tail=[tail_extension]
+```
+
+If the equilibration estimator is used, `_eq` is appended to the filename.

scripts/poisson.jl

@@ -4,7 +4,7 @@
 
 runs AFEM loop for stochastic Poisson problem
  
-main usage:
+usage:
 - run experiment: run(; problem = problem, kwargs...)
 - load results  : show_results(; kwargs...)
 - produce plots : produce_plots(; kwargs...)
@@ -78,6 +78,42 @@ default_args = Dict(
 )
 
 
+## that is the main function to run and save results
+function run(; force = false, kwargs...)
+    data = deepcopy(default_args)
+
+    for (k, v) in kwargs
+        data[String(k)] = v
+    end
+    Plotter = data["Plotter"]
+
+    ## load/produce data
+    data, ~ = produce_or_load(_main, data, filename = filename, force = force)
+
+    ## plot final solution
+    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
+
+
+# loads (or produces) a result file and prints the data
 function show_results(; force = false, mode_history_up_to_level = 15, kwargs...)
     data = deepcopy(default_args)
 
@@ -113,8 +149,9 @@ function show_results(; force = false, mode_history_up_to_level = 15, kwargs...)
     return data
 end
 
-
-function main(
+# the main AFEM loop for the given problem and data
+# (this function does not check for existing results, better use run !)
+function _main(
         data = nothing;
         debug = false,
         plot_solution = false,
@@ -396,40 +433,6 @@ function main(
 end
 
 
-function run(; force = false, kwargs...)
-    data = deepcopy(default_args)
-
-    for (k, v) in kwargs
-        data[String(k)] = v
-    end
-    Plotter = data["Plotter"]
-
-    ## load/produce data
-    data, ~ = produce_or_load(main, data, filename = filename, force = force)
-
-    ## plot final solution
-    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 produce_plots(;
         order = default_args["order"],
         decay = default_args["decay"],

src/modelproblems/logpoisson_dual.jl

@@ -1,7 +1,24 @@
 """
 $(TYPEDEF)
 
-Dual log-transformed formulation of the Poisson problem with exponential stochastic coefficient (WIP).
+
+Poisson problem with exponential stochastic coefficient ``e^a`` that seeks ``u`` such that
+
+``-\\mathrm{div}(e^a(y,x) \\nabla u(y,x)) = f(x) \\quad \\text{for } (y,x) \\in \\Gamma \\times D``
+
+The dual formulation of the log-transformed formulation of the Poisson problem
+introduces the auxiliary stress variable ``p := - e^a(y,x) \\nabla u(y,x) = -∇ũ - ∇aũ`` for the transformed
+``ũ := e^{-a} u``.
+
+Its weak formulation seeks ``(p,ũ)`` such that
+
+```math
+\\begin{aligned}
+(p, q) + (ũ \\nabla a, q) - (\\mathrm{div}(q), ũ) & = 0 \\quad \\text{for all } q \\\\
+-(\\mathrm{div}(p), v) & = (f, v) \\quad \\text{for all } v
+\\end{aligned}
+```
+
 """
 abstract type LogTransformedPoissonProblemDual <: AbstractModelProblem end