From 1d24f154bb3de2ddc4531848076c6c6a57ef4cd9 Mon Sep 17 00:00:00 2001 From: xtalax Date: Thu, 2 Mar 2023 16:29:27 +0000 Subject: [PATCH 01/26] functionscheme, params -> ivs --- docs/pages.jl | 2 +- docs/src/api/utils.md | 2 +- docs/src/devnotes.md | 2 +- docs/src/tutorials/params.md | 4 +- src/MOL_discretization.jl | 2 +- src/MOL_utils.jl | 2 +- .../array_form/stencil_utils.jl | 174 ++++++++++++++++++ src/discretization/discretize_vars.jl | 22 ++- src/discretization/generate_bc_eqs.jl | 8 +- src/discretization/generate_ic_defaults.jl | 2 +- src/discretization/schemes/WENO/WENO.jl | 2 +- src/discretization/schemes/WENO/WENO_array.jl | 89 +++++++++ .../centered_difference.jl | 2 +- .../function_scheme/function_scheme.jl | 38 ++++ .../integral_expansion/integral_expansion.jl | 6 +- .../integral_expansion_array.jl | 46 +++++ .../nonlinear_laplacian.jl | 12 +- .../spherical_laplacian.jl | 6 +- .../upwind_difference/upwind_difference.jl | 6 +- .../upwind_difference_array.jl | 135 ++++++++++++++ src/interface/scheme_types.jl | 49 +++++ src/system_parsing/bcs/parse_boundaries.jl | 4 +- src/system_parsing/variable_map.jl | 16 +- test/components/finite_diff_schemes.jl | 32 ++-- test/components/function_scheme.jl | 11 ++ 25 files changed, 611 insertions(+), 63 deletions(-) create mode 100644 src/discretization/array_form/stencil_utils.jl create mode 100644 src/discretization/schemes/WENO/WENO_array.jl create mode 100644 src/discretization/schemes/function_scheme/function_scheme.jl create mode 100644 src/discretization/schemes/integral_expansion/integral_expansion_array.jl create mode 100644 src/discretization/schemes/upwind_difference/upwind_difference_array.jl create mode 100644 test/components/function_scheme.jl diff --git a/docs/pages.jl b/docs/pages.jl index a7c55f7f1..32c17501e 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -3,7 +3,7 @@ pages = [ "Tutorials" => [ "tutorials/brusselator.md", "tutorials/heat.md", - "tutorials/params.md", + "tutorials/ivs.md", "tutorials/heatss.md", "tutorials/sispde.md", "tutorials/icbc_sampled.md", diff --git a/docs/src/api/utils.md b/docs/src/api/utils.md index d8d944dac..8c886f111 100644 --- a/docs/src/api/utils.md +++ b/docs/src/api/utils.md @@ -6,6 +6,6 @@ CurrentModule = MethodOfLines ```@docs unitindex -params +ivs Idx ``` \ No newline at end of file diff --git a/docs/src/devnotes.md b/docs/src/devnotes.md index 54f252420..7bb5698f0 100644 --- a/docs/src/devnotes.md +++ b/docs/src/devnotes.md @@ -48,7 +48,7 @@ end # Note that indexmap is used along with the function `Idx` to create an equivalent index for the discrete form of `u`, # which may have a different number of dimensions to `II` function generate_central_difference_rules(II::CartesianIndex, s::DiscreteSpace, terms::Vector{<:Term}, indexmap::Dict) - rules = [[@rule Differential(x)(u) => second_order_central_difference(Idx(II, s, u, indexmap), s, u, x) for x in params(u, s)] for u in depvars] + rules = [[@rule Differential(x)(u) => second_order_central_difference(Idx(II, s, u, indexmap), s, u, x) for x in ivs(u, s)] for u in depvars] rules = reduce(vcat, rules) diff --git a/docs/src/tutorials/params.md b/docs/src/tutorials/params.md index bd044c243..882c375e0 100644 --- a/docs/src/tutorials/params.md +++ b/docs/src/tutorials/params.md @@ -2,7 +2,7 @@ We can also build up more complicated systems with multiple dependent variables and parameters as follows -```@example params1 +```@example ivs1 using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets @parameters t x @@ -50,7 +50,7 @@ gif(anim, "plot.gif",fps=30) The system does not need to be re-discretized every time we want to plot with different parameters, the system can be remade with new parameters with `remake`. See the `ModelingToolkit.jl` [docs](https://docs.sciml.ai/ModelingToolkit/stable/tutorials/ode_modeling/#Algebraic-relations-and-structural-simplification) for more ways to manipulate a `prob` post discretization. -```@example params2 +```@example ivs2 using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets @parameters t x diff --git a/src/MOL_discretization.jl b/src/MOL_discretization.jl index 38f0ccab1..e884ef99e 100644 --- a/src/MOL_discretization.jl +++ b/src/MOL_discretization.jl @@ -107,7 +107,7 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method eqvar = interiormap.var[pde] # * Assumes that all variables in the equation have same dimensionality except edgevals - args = params(eqvar, s) + args = ivs(eqvar, s) indexmap = Dict([args[i] => i for i in 1:length(args)]) if disc_strategy isa ScalarizedDiscretization # Generate the equations for the interior points diff --git a/src/MOL_utils.jl b/src/MOL_utils.jl index c81418190..9cb38d0bd 100644 --- a/src/MOL_utils.jl +++ b/src/MOL_utils.jl @@ -100,7 +100,7 @@ end function newindex(u_, II, s, indexmap) u = depvar(u_, s) args_ = remove(arguments(u_), s.time) - args = params(u, s) + args = ivs(u, s) is = map(enumerate(args_)) do (j, x) if haskey(indexmap, x) II[indexmap[x]] diff --git a/src/discretization/array_form/stencil_utils.jl b/src/discretization/array_form/stencil_utils.jl new file mode 100644 index 000000000..7009dcce2 --- /dev/null +++ b/src/discretization/array_form/stencil_utils.jl @@ -0,0 +1,174 @@ +function lower_boundary_deriv(D, udisc, iboundary, j, is, interior) + weights = D.low_boundary_coefs[iboundary] + taps = 1:D.boundary_stencil_length + prepare_boundary_op((BoundaryDerivArrayOp(weights, taps, udisc, j, is, interior), + iboundary), interior, j) +end + +function upper_boundary_deriv(D, udisc, iboundary, j, is, interior, lenx) + weights = D.high_boundary_coefs[lenx - iboundary + 1] + taps = (lenx-D.boundary_stencil_length+1):lenx + prepare_boundary_op((BoundaryDerivArrayOp(weights, taps, udisc, j, is, interior), + iboundary), interior, j) +end + +function integral_op_pair(dx, udisc, j, is, ranges, interior, i) + prepare_boundary_op((IntegralArrayOp(dx, udisc, i, j, is, ranges, interior), + i), interior, j) +end + +function prepare_boundary_op(boundaryop, interior, j) + function maketuple(i) + out = map(1:length(interior)) do k + k == j ? i : interior[k] + end + return Tuple(out) + end + (op, iboundary) = boundaryop + return maketuple(iboundary) => op +end + +function prepare_boundary_ops(boundaryops, interior, j) + function maketuple(i) + out = map(1:length(interior)) do k + k == j ? i : interior[k] + end + return Tuple(out) + end + return map(boundaryops) do (op, iboundary) + maketuple(iboundary) => op + end +end + +function interior_deriv(D::DerivativeOperator{T,N,Wind,DX}, udisc, s, offsets, j, is, interior, bs, isx = false) where {T,N,Wind,DX<:Number} + weights = D.stencil_coefs + taps = offsets .+ is[j] + InteriorDerivArrayOp(weights, taps, udisc, s, j, is, interior, bs, isx) +end + +function interior_deriv(D::DerivativeOperator{T,N,Wind,DX}, udisc, s, offsets, j, is, interior, bs, isx = false) where {T,N,Wind,DX<:AbstractVector} + @assert length(bs) == 0 "Interface boundary conditions are not yet supported for nonuniform dx dimensions, such as $x, please post an issue to https://github.com/SciML/MethodOfLines.jl if you need this functionality." + weights = D.stencil_coefs[is[j]-D.boundary_point_count] + taps = offsets .+ is[j] + InteriorDerivArrayOp(weights, taps, udisc, s, j, is, interior, bs, isx) +end + +function BoundaryDerivArrayOp(weights, taps, udisc, j, is, interior) + # * I Possibly needs updating + Is = map(taps) do tap + map(1:ndims(udisc)) do i + if i == j + tap + else + is[i] + end + end + end + expr = sym_dot(weights, map(I -> udisc[I...], Is)) + + symindices = setdiff(1:ndims(udisc), [j]) + output_idx = Tuple(is[symindices]) + return FillArrayOp(recursive_unwrap(expr), output_idx, interior[symindices]) +end + +reduce_interior(interior, j) = first(interior[j]) == 1 ? [interior[1:j-1]..., 2:last(interior[j]), interior[j+1:end]...] : interior + +function trapezium_sum(ranges, interior, dx, udisc, is, j) + N = ndims(udisc) + I1 = unitindex(N, j) + I = CartesianIndex(is...) + i = is[j] + Im1 = I - I1 + im1 = is[j] - 1 + rinterior = reduce_interior(interior, j) + if dx isa Number + expr = (dx*(udisc[Im1] + udisc[I]) / 2) + else + expr = (dx[im1]*udisc[Im1] + dx[i]*udisc[I]) / 2 + end + + return FillArrayMaker(expr, is, ranges, rinterior)[interior...] +end + +function IntegralArrayOp(dx, udisc, k, j, is, ranges, interior, iswd = false) + if iswd + interior = [interior[1:j-1]..., 1:size(udisc, j), interior[j+1:end]...] + end + trapop = trapezium_sum(ranges, interior, dx, udisc, is, j) + op = sum(trapop[1:(k-first(interior[j])+1)], dims = j) + + return op +end + +function IntegralArrayMaker(dx, udisc, k, j, is, ranges, interior, iswd = false, ) + if iswd + ranges = [ranges[1:j-1]..., 1:size(udisc, j), ranges[j+1:end]...] + end + + op = IntegralArrayOp(dx, udisc, k, j, is, ranges, interior, iswd) + return FillArrayMaker(op, is, ranges, interior) +end + +function InteriorDerivArrayOp(weights, taps, udisc, s, j, output_idx, interior, bs, isx = false) + # * I Possibly needs updating + Is = map(taps) do tap + _is = map(1:ndims(udisc)) do i + if i == j + tap + else + output_idx[i] + end + end + CartesianIndex(_is...) + end + # Wrap interfaces + Is = map(I -> bwrap(I, bs, s, j, isx), Is) + + Is = map(Is) do I + map(1:ndims(udisc)) do i + I[i] + end + end + + expr = sym_dot(weights, map(I -> udisc[I...], Is)) + + return FillArrayOp(expr, output_idx, interior) +end + +function FillArrayOp(expr, output_idx, interior) + ranges = Dict(output_idx .=> interior) # hope this doesn't check bounds eagerly + return ArrayOp(Array{symtype(expr),length(output_idx)}, + output_idx, expr, +, nothing, ranges) +end + +NullBG_ArrayMaker(ranges, ops) = ArrayMaker{Real}(Tuple(map(r -> r[end] - r[1] + 1, ranges)), + vcat(Tuple(ranges) => 0, ops)) +Construct_ArrayMaker(ranges, + ops) = ArrayMaker{Real}(Tuple(map(r -> r[end] - r[1] + 1, ranges)), ops) +Construct_ArrayMaker{T}(ranges, + ops) = ArrayMaker{T}(Tuple(map(r -> r[end] - r[1] + 1, ranges)), ops) + + + FillArrayMaker(expr, is, + ranges, interior) = NullBG_ArrayMaker(ranges, + [Tuple(interior) => FillArrayOp(expr, is, interior)]) + +ArrayMakerWrap(udisc, ranges) = Arraymaker{Real}(Tuple(map(r -> r[end] - r[1] + 1, ranges)), + [Tuple(ranges) => udisc]) + +##### + +function get_interior(u, s, interior) + map(ivs(u, s)) do x + if haskey(interior, x) + interior[x] + else + 1:length(s, x) + end + end +end + +function get_ranges(u, s) + map(x -> first(axes(s.grid[x])), ivs(u, s)) +end +get_is(u, s) = map(x -> s.index_syms[x], ivs(u, s)) diff --git a/src/discretization/discretize_vars.jl b/src/discretization/discretize_vars.jl index 6f75312d3..271754df4 100644 --- a/src/discretization/discretize_vars.jl +++ b/src/discretization/discretize_vars.jl @@ -149,28 +149,30 @@ function DiscreteSpace(vars, discretization::MOLFiniteDifference{G}) where {G} end function Base.getproperty(s::DiscreteSpace, p::Symbol) - if p in [:ū, :x̄, :time, :args, :x2i, :i2x] + if p in [:ū, :x̄, :ps, :time, :args, :x2i, :i2x] getfield(s.vars, p) else getfield(s, p) end end +params(s::DiscreteSpace) = s.ps + get_grid_type(::DiscreteSpace{N,M,G}) where {N,M,G} = G prepare_dx(dx::Integer, xdomain, ::CenterAlignedGrid) = (xdomain[2] - xdomain[1])/(dx - 1) prepare_dx(dx::Integer, xdomain, ::EdgeAlignedGrid) = (xdomain[2] - xdomain[1])/dx prepare_dx(dx, xdomain, ::AbstractGrid) = dx -nparams(::DiscreteSpace{N,M}) where {N,M} = N +nivs(::DiscreteSpace{N,M}) where {N,M} = N nvars(::DiscreteSpace{N,M}) where {N,M} = M """ - params(u, s::DiscreteSpace) + ivs(u, s::DiscreteSpace) Filter out the time variable and get the spatial variables of `u` in `s`. """ -params(u, s::DiscreteSpace) = remove(s.args[operation(u)], s.time) +ivs(u, s::DiscreteSpace) = remove(s.args[operation(u)], s.time) Base.ndims(u, s::DiscreteSpace) = ndims(s.discvars[depvar(u, s)]) Base.length(s::DiscreteSpace, x) = length(s.grid[x]) @@ -185,9 +187,9 @@ of `II` that corresponds to only the spatial arguments of `u`. """ function Idx(II::CartesianIndex, s::DiscreteSpace, u, indexmap) # We need to construct a new index as indices may be of different size - length(params(u, s)) == 0 && return CartesianIndex() - !all(x -> haskey(indexmap, x), params(u, s)) && return II - is = [II[indexmap[x]] for x in params(u, s)] + length(ivs(u, s)) == 0 && return CartesianIndex() + !all(x -> haskey(indexmap, x), ivs(u, s)) && return II + is = [II[indexmap[x]] for x in ivs(u, s)] II = CartesianIndex(is...) return II @@ -198,7 +200,7 @@ A function that returns what to replace independent variables with in boundary e """ @inline function axiesvals(s::DiscreteSpace{N,M,G}, u_, x_, I) where {N,M,G} u = depvar(u_, s) - map(params(u, s)) do x + map(ivs(u, s)) do x if isequal(x, x_) x => (I[x2i(s, u, x)] == 1 ? first(s.axies[x]) : last(s.axies[x])) else @@ -207,8 +209,8 @@ A function that returns what to replace independent variables with in boundary e end end -gridvals(s::DiscreteSpace{N}, u) where {N} = ndims(u, s) == 0 ? [] : map(y -> [x => s.grid[x][y.I[x2i(s, u, x)]] for x in params(u, s)], s.Igrid[u]) -gridvals(s::DiscreteSpace{N}, u, I::CartesianIndex) where {N} = ndims(u, s) == 0 ? [] : [x => s.grid[x][I[x2i(s, u, x)]] for x in params(u, s)] +gridvals(s::DiscreteSpace{N}, u) where {N} = ndims(u, s) == 0 ? [] : map(y -> [x => s.grid[x][y.I[x2i(s, u, x)]] for x in ivs(u, s)], s.Igrid[u]) +gridvals(s::DiscreteSpace{N}, u, I::CartesianIndex) where {N} = ndims(u, s) == 0 ? [] : [x => s.grid[x][I[x2i(s, u, x)]] for x in ivs(u, s)] varmaps(s::DiscreteSpace, depvars, II, indexmap) = [u => s.discvars[u][Idx(II, s, u, indexmap)] for u in depvars] diff --git a/src/discretization/generate_bc_eqs.jl b/src/discretization/generate_bc_eqs.jl index 769826c03..f94d508ef 100644 --- a/src/discretization/generate_bc_eqs.jl +++ b/src/discretization/generate_bc_eqs.jl @@ -1,5 +1,5 @@ @inline function generate_bc_eqs!(bceqs, s, boundaryvalfuncs, interiormap, boundary::AbstractTruncatingBoundary) - args = params(depvar(boundary.u, s), s) + args = ivs(depvar(boundary.u, s), s) indexmap = Dict([args[i]=>i for i in 1:length(args)]) push!(bceqs, generate_bc_eqs(s, boundaryvalfuncs, boundary, interiormap, indexmap)) end @@ -48,7 +48,7 @@ function boundary_value_maps(II, s::DiscreteSpace{N,M,G}, boundary, derivweights # replace u(t,0) with u₁, etc u = depvar(u_, s) - args = params(u, s) + args = ivs(u, s) j = findfirst(isequal(x_), args) IIold = II # We need to construct a new index in case the value at the boundary appears in an equation one dimension lower @@ -99,7 +99,7 @@ function boundary_value_maps(II, s::DiscreteSpace{N,M,G}, boundary, derivweights # * Assume that the BC is in terms of an explicit expression, not containing references to variables other than u_ at the boundary u = depvar(u_, s) - args = params(u, s) + args = ivs(u, s) j = findfirst(isequal(x_), args) IIold = II # We need to construct a new index in case the value at the boundary appears in an equation one dimension lower @@ -155,7 +155,7 @@ Pads the boundaries with extrapolation equations, extrapolated with 6th order la Reuses `central_difference` as this already dispatches the correct stencil, given a `DerivativeOperator` which contains the correct weights. """ function generate_extrap_eqs!(eqs, pde, u, s, derivweights, interiormap, bcmap) - args = params(u, s) + args = ivs(u, s) length(args) == 0 && return lowerextents, upperextents = interiormap.stencil_extents[pde] diff --git a/src/discretization/generate_ic_defaults.jl b/src/discretization/generate_ic_defaults.jl index 8b0e10fe6..d922bfa62 100644 --- a/src/discretization/generate_ic_defaults.jl +++ b/src/discretization/generate_ic_defaults.jl @@ -6,7 +6,7 @@ function generate_ic_defaults(tconds, s, ::ScalarizedDiscretization) throw(ArgumentError("Upper boundary condition $(ic.eq) on time variable is not supported, please use a change of variables `t => -τ` to make this an initial condition.")) end - args = params(depvar(ic.u, s), s) + args = ivs(depvar(ic.u, s), s) indexmap = Dict([args[i]=>i for i in 1:length(args)]) D = ic.order == 0 ? identity : (Differential(t)^ic.order) defaultvars = D.(s.discvars[depvar(ic.u, s)]) diff --git a/src/discretization/schemes/WENO/WENO.jl b/src/discretization/schemes/WENO/WENO.jl index 17ba3fd9c..827d0e72f 100644 --- a/src/discretization/schemes/WENO/WENO.jl +++ b/src/discretization/schemes/WENO/WENO.jl @@ -85,5 +85,5 @@ end This is a catch all ruleset, as such it does not use @rule. """ @inline function generate_WENO_rules(II::CartesianIndex, s::DiscreteSpace, depvars, derivweights::DifferentialDiscretizer, bcmap, indexmap, terms) - return reduce(safe_vcat, [[(Differential(x))(u) => weno(Idx(II, s, u, indexmap), s, derivweights.advection_scheme, filter_interfaces(bcmap[operation(u)][x]), (x2i(s, u, x), x), u, s.dxs[x]) for x in params(u, s)] for u in depvars], init = []) + return reduce(safe_vcat, [[(Differential(x))(u) => weno(Idx(II, s, u, indexmap), s, derivweights.advection_scheme, filter_interfaces(bcmap[operation(u)][x]), (x2i(s, u, x), x), u, s.dxs[x]) for x in ivs(u, s)] for u in depvars], init = []) end diff --git a/src/discretization/schemes/WENO/WENO_array.jl b/src/discretization/schemes/WENO/WENO_array.jl new file mode 100644 index 000000000..31012fcc4 --- /dev/null +++ b/src/discretization/schemes/WENO/WENO_array.jl @@ -0,0 +1,89 @@ +######################################################################################## +# Stencil interface +######################################################################################## + +function weno(interior, s::DiscreteSpace, wenoscheme::WENOScheme, bs, jx, u, dx::Number) + j, x = jx + ε = wenoscheme.epsilon + + interior = get_interior(u, s, interior) + ranges = get_ranges(u, s) + is = get_is(u, s) + + II = CartesianIndex(is...) + I1 = unitindex(ndims(u, s), j) + + udisc = s.discvars[u] + + Im2 = bwrap(II - 2I1, bs, s, j) + Im1 = bwrap(II - I1, bs, s, j) + Ip1 = bwrap(II + I1, bs, s, j) + Ip2 = bwrap(II + 2I1, bs, s, j) + + u_m2 = udisc[Im2] + u_m1 = udisc[Im1] + u_0 = udisc[II] + u_p1 = udisc[Ip1] + u_p2 = udisc[Ip2] + + γm1 = 1 / 10 + γm2 = 3 / 5 + γm3 = 3 / 10 + + β1 = 13 * (u_0 - 2 * u_p1 + u_p2)^2 / 12 + (3 * u_0 - 4 * u_p1 + u_p2)^2 / 4 + β2 = 13 * (u_m1 - 2 * u_0 + u_p1)^2 / 12 + (u_m1 - u_p1)^2 / 4 + β3 = 13 * (u_m2 - 2 * u_m1 + u_0)^2 / 12 + (u_m2 - 4 * u_m1 + 3 * u_0)^2 / 4 + + ωm1 = γm1 / (ε + β1)^2 + ωm2 = γm2 / (ε + β2)^2 + ωm3 = γm3 / (ε + β3)^2 + + wm_denom = ωm1 + ωm2 + ωm3 + wm1 = ωm1 / wm_denom + wm2 = ωm2 / wm_denom + wm3 = ωm3 / wm_denom + + γp1 = 3 / 10 + γp2 = 3 / 5 + γp3 = 1 / 10 + + ωp1 = γp1 / (ε + β1)^2 + ωp2 = γp2 / (ε + β2)^2 + ωp3 = γp3 / (ε + β3)^2 + + wp_denom = ωp1 + ωp2 + ωp3 + wp1 = ωp1 / wp_denom + wp2 = ωp2 / wp_denom + wp3 = ωp3 / wp_denom + + hm1 = (11 * u_0 - 7 * u_p1 + 2 * u_p2) / 6 + hm2 = (5 * u_0 - u_p1 + 2 * u_m1) / 6 + hm3 = (2 * u_0 + 5 * u_m1 - u_m2) / 6 + + hp1 = (2 * u_0 + 5 * u_p1 - u_p2) / 6 + hp2 = (5 * u_0 + 2 * u_p1 - u_m1) / 6 + hp3 = (11 * u_0 - 7 * u_m1 + 2 * u_m2) / 6 + + hp = wp1 * hp1 + wp2 * hp2 + wp3 * hp3 + hm = wm1 * hm1 + wm2 * hm2 + wm3 * hm3 + + expr = (hp - hm) / dx + return FillArrayMaker(recursive_unwrap(expr), is, ranges, interior)[interior...] +end + +function weno(interior, s::DiscreteSpace, b, jx, u, dx::AbstractVector) + @assert false "WENO scheme not implemented for nonuniform grids." +end + +""" +This is a catch all ruleset, as such it does not use @rule. +""" +@inline function generate_WENO_rules(interior, s::DiscreteSpace, depvars, + derivweights::DifferentialDiscretizer, pmap, + indexmap, terms) + return reduce(safe_vcat, + [[(Differential(x))(u) => weno(interior, s, derivweights.advection_scheme, + pmap.map[operation(u)][x], (x2i(s, u, x), x), u, s.dxs[x]) + for x in ivs(u, s)] + for u in depvars], init = []) +end diff --git a/src/discretization/schemes/centered_difference/centered_difference.jl b/src/discretization/schemes/centered_difference/centered_difference.jl index 05aab1dc9..5869e88b4 100644 --- a/src/discretization/schemes/centered_difference/centered_difference.jl +++ b/src/discretization/schemes/centered_difference/centered_difference.jl @@ -60,5 +60,5 @@ This is a catch all ruleset, as such it does not use @rule. Any even ordered der let orders = derivweights.orders[x] orders[iseven.(orders)] end - )] for x in params(u, s)], init = []) for u in depvars], init = []) + )] for x in ivs(u, s)], init = []) for u in depvars], init = []) end diff --git a/src/discretization/schemes/function_scheme/function_scheme.jl b/src/discretization/schemes/function_scheme/function_scheme.jl new file mode 100644 index 000000000..b56078c1f --- /dev/null +++ b/src/discretization/schemes/function_scheme/function_scheme.jl @@ -0,0 +1,38 @@ +function function_scheme(F::FunctionScheme{is_nonuniform}, II, s, bs, jx, u, ufunc) where {is_nonuniform <: Val{false}} + + + # Tap points of the stencil, this uses boundary_point_count as this is equal to half the stencil size, which is what we want. + u_disc = ufunc(u, Itap, x) + ps = params(s) + dx = s.dxs[x] + t = s.time + + return f(u_disc, ps, t, x, dx) +end + + +function get_f_and_taps(F::FunctionScheme{is_nonuniform}, II, s, bs, jx, u, ufunc) where {is_nonuniform <: Val{false}} + j, x = jx + ndims(u, s) == 0 && return 0 + # unit index in direction of the derivative + I1 = unitindex(ndims(u, s), j) + # offset is important due to boundary proximity + haslower, hasupper = haslowerupper(bs, x) + + lower_point_count = length(f.lower) + upper_point_count = length(f.upper) + + if (II[j] <= D.lower_point_count) & !haslower + f = F.lower[II[j]] + offset = 1 - II[j] + Itap = [II + (i + offset) * I1 for i in 0:(F.boundary_points-1)] + elseif (II[j] > (length(s, x) - upper_point_count)) & !hasupper + f = F.upper[length(s, x)-II[j]+1] + offset = length(s, x) - II[j] + Itap = [II + (i + offset) * I1 for i in (-F.boundary_points+1):1:0] + else + f = F.interior + Itap = [bwrap(II + i * I1, bs, s, jx) for i in half_range(F.interior_points)] + end + + return f, Itap diff --git a/src/discretization/schemes/integral_expansion/integral_expansion.jl b/src/discretization/schemes/integral_expansion/integral_expansion.jl index 8b0e5c82d..f7bccf73b 100644 --- a/src/discretization/schemes/integral_expansion/integral_expansion.jl +++ b/src/discretization/schemes/integral_expansion/integral_expansion.jl @@ -52,14 +52,14 @@ end ufunc(u, I, x) = s.discvars[u][I] eulerrules = reduce(safe_vcat, [[Integral(x in DomainSets.ClosedInterval(s.vars.intervals[x][1], Num(x)))(u) => euler_integral(Idx(II, s, u, indexmap), s, (x2i(s, u, x), x), u, ufunc) - for x in params(u, s)] + for x in ivs(u, s)] for u in depvars], init = []) return eulerrules end function wd_integral_Idx(II::CartesianIndex, s::DiscreteSpace, u, x, indexmap) # We need to construct a new index as indices may be of different size - length(params(u, s)) == 0 && return CartesianIndex() + length(ivs(u, s)) == 0 && return CartesianIndex() # A hack using the boundary value re-indexing function to get an index that will work u_ = substitute(u, [x => s.axies[x][end]]) II = newindex(u_, II, s, indexmap) @@ -71,7 +71,7 @@ end wholedomainrules = reduce(safe_vcat, [[Integral(x in DomainSets.ClosedInterval(s.vars.intervals[x][1], s.vars.intervals[x][2]))(u) => whole_domain_integral(wd_integral_Idx(II, s, u, x, indexmap), s, (x2i(s, u, x), x), u, ufunc) - for x in filter(x -> (!haskey(indexmap, x) | isequal(x, bvar)), params(u, s))] + for x in filter(x -> (!haskey(indexmap, x) | isequal(x, bvar)), ivs(u, s))] for u in depvars], init = []) return wholedomainrules diff --git a/src/discretization/schemes/integral_expansion/integral_expansion_array.jl b/src/discretization/schemes/integral_expansion/integral_expansion_array.jl new file mode 100644 index 000000000..d10f5b034 --- /dev/null +++ b/src/discretization/schemes/integral_expansion/integral_expansion_array.jl @@ -0,0 +1,46 @@ +function euler_integral(interior, s, jx, u, udisc) + j, x = jx + dx = s.dxs[x] + interior = get_interior(u, s, interior) + ranges = get_ranges(u, s) + is = get_is(u, s) + + oppairs = map(interior[j]) do i + integral_op_pair(dx, udisc, j, is, ranges, interior, i) + end + + return NullBG_ArrayMaker(ranges, oppairs)[interior...] +end + +# An integral across the whole domain (xmin .. xmax) +function whole_domain_integral(interior, s, jx, u, udisc) + j, x = jx + dx = s.dxs[x] + + interior = get_interior(u, s, interior) + ranges = get_ranges(u, s) + is = get_is(u, s) + + lenx = length(s, x) + + return IntegralArrayMaker(dx, udisc, lenx, j, is, ranges, interior, true)[interior...] +end + +@inline function generate_euler_integration_rules(interior, s::DiscreteSpace, depvars, indexmap, terms) + eulerrules = reduce(safe_vcat, [[Integral(x in DomainSets.ClosedInterval(s.vars.intervals[x][1], Num(x)))(u) => + euler_integral(interior, s, (x2i(s, u, x), x), u, s.discvars[u]) + for x in ivs(u, s)] + for u in depvars], init = []) + return eulerrules +end + +@inline function generate_whole_domain_integration_rules(interior, s::DiscreteSpace, depvars, indexmap, terms, bvar = nothing) + wholedomainrules = reduce(safe_vcat, + [[Integral(x in DomainSets.ClosedInterval(s.vars.intervals[x][1], s.vars.intervals[x][2]))(u) => + whole_domain_integral(interior, s, (x2i(s, u, x), x), u, s.discvars[u]) + for x in filter(x -> (!haskey(indexmap, x) | isequal(x, bvar)), ivs(u, s))] + for u in depvars], + init = []) + + return wholedomainrules +end diff --git a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl index d97ca7886..78d0db9dd 100644 --- a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl +++ b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl @@ -55,7 +55,7 @@ function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace, # Map parameters to interpolated values. Using simplistic extrapolation/interpolation for now as grids are uniform #TODO: make this more efficient - map_params_to_interpolated(stencil, weights) = safe_vcat([x => dot(weights, getindex.((s.grid[x],), getindex.(interface_wrap(stencil), (j,))))], [s.x̄[k] => s.grid[s.x̄[k]][II[k]] for k in setdiff(1:N, [j])]) + map_ivs_to_interpolated(stencil, weights) = safe_vcat([x => dot(weights, getindex.((s.grid[x],), getindex.(interface_wrap(stencil), (j,))))], [s.x̄[k] => s.grid[s.x̄[k]][II[k]] for k in setdiff(1:N, [j])]) # Take the inner finite difference inner_difference = [dot(inner_weights, s.discvars[u][interface_wrap(inner_stencil)]) for (inner_weights, inner_stencil) in inner_deriv_weights_and_stencil] @@ -64,7 +64,7 @@ function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace, interpolated_expr = map(interp_weights_and_stencil) do (weights, stencil) - substitute(substitute(expr, map_vars_to_interpolated(stencil, weights)), map_params_to_interpolated(stencil, weights)) + substitute(substitute(expr, map_vars_to_interpolated(stencil, weights)), map_ivs_to_interpolated(stencil, weights)) end # multiply the inner finite difference by the interpolated expression, and finally take the outer finite difference @@ -72,13 +72,13 @@ function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace, end @inline function generate_nonlinlap_rules(II::CartesianIndex, s::DiscreteSpace, depvars, derivweights::DifferentialDiscretizer, bcmap, indexmap, terms) - rules = reduce(safe_vcat, [vec([@rule *(~~c, $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)), ~~d) => *(~c..., cartesian_nonlinear_laplacian(*(~a..., ~b...), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u), ~d...) for x in params(u, s)]) for u in depvars], init = []) + rules = reduce(safe_vcat, [vec([@rule *(~~c, $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)), ~~d) => *(~c..., cartesian_nonlinear_laplacian(*(~a..., ~b...), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u), ~d...) for x in ivs(u, s)]) for u in depvars], init = []) - rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)) => cartesian_nonlinear_laplacian(*(~a..., ~b...), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u) for x in params(u, s)]) for u in depvars], init = [])) + rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)) => cartesian_nonlinear_laplacian(*(~a..., ~b...), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u) for x in ivs(u, s)]) for u in depvars], init = [])) - rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule ($(Differential(x))($(Differential(x))(u) / ~a)) => cartesian_nonlinear_laplacian(1 / ~a, Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u) for x in params(u, s)]) for u in depvars], init = [])) + rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule ($(Differential(x))($(Differential(x))(u) / ~a)) => cartesian_nonlinear_laplacian(1 / ~a, Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u) for x in ivs(u, s)]) for u in depvars], init = [])) - rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule *(~~b, ($(Differential(x))($(Differential(x))(u) / ~a)), ~~c) => *(~b..., ~c..., cartesian_nonlinear_laplacian(1 / ~a, Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u)) for x in params(u, s)]) for u in depvars], init = [])) + rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule *(~~b, ($(Differential(x))($(Differential(x))(u) / ~a)), ~~c) => *(~b..., ~c..., cartesian_nonlinear_laplacian(1 / ~a, Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u)) for x in ivs(u, s)]) for u in depvars], init = [])) nonlinlap_rules = [] for t in terms diff --git a/src/discretization/schemes/spherical_laplacian/spherical_laplacian.jl b/src/discretization/schemes/spherical_laplacian/spherical_laplacian.jl index 013a752e6..9fcb666f1 100644 --- a/src/discretization/schemes/spherical_laplacian/spherical_laplacian.jl +++ b/src/discretization/schemes/spherical_laplacian/spherical_laplacian.jl @@ -36,13 +36,13 @@ end @inline function generate_spherical_diffusion_rules(II::CartesianIndex, s::DiscreteSpace, depvars, derivweights::DifferentialDiscretizer, bcmap, indexmap, terms) rules = reduce(safe_vcat, [vec([@rule *(~~a, 1 / (r^2), ($(Differential(r))(*(~~c, (r^2), ~~d, $(Differential(r))(u), ~~e))), ~~b) => *(~a..., spherical_diffusion(*(~c..., ~d..., ~e..., Num(1)), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][r]), depvars, r, u), ~b...) - for r in params(u, s)]) for u in depvars], init = []) + for r in ivs(u, s)]) for u in depvars], init = []) rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule /(*(~~a, $(Differential(r))(*(~~c, (r^2), ~~d, $(Differential(r))(u), ~~e)), ~~b), (r^2)) => *(~a..., ~b..., spherical_diffusion(*(~c..., ~d..., ~e..., Num(1)), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][r]), depvars, r, u)) - for r in params(u, s)]) for u in depvars], init = [])) + for r in ivs(u, s)]) for u in depvars], init = [])) rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule /(($(Differential(r))(*(~~c, (r^2), ~~d, $(Differential(r))(u), ~~e))), (r^2)) => spherical_diffusion(*(~c..., ~d..., ~e..., Num(1)), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][r]), depvars, r, u) - for r in params(u, s)]) for u in depvars], init = [])) + for r in ivs(u, s)]) for u in depvars], init = [])) spherical_diffusion_rules = [] for t in terms diff --git a/src/discretization/schemes/upwind_difference/upwind_difference.jl b/src/discretization/schemes/upwind_difference/upwind_difference.jl index 1feab1bdd..60071e5cb 100644 --- a/src/discretization/schemes/upwind_difference/upwind_difference.jl +++ b/src/discretization/schemes/upwind_difference/upwind_difference.jl @@ -84,14 +84,14 @@ end let orders = derivweights.orders[x] setdiff(orders[isodd.(orders)], skip) end - )] for x in params(u, s)], init = []) for u in depvars], init = []), + )] for x in ivs(u, s)], init = []) for u in depvars], init = []), #Catch division and multiplication, see issue #1 reduce(safe_vcat, [reduce(safe_vcat, [[@rule /(*(~~a, $(Differential(x)^d)(u), ~~b), ~c) => upwind_difference(*(~a..., ~b...) / ~c, d, Idx(II, s, u, indexmap), s, filter_interfaces(bcmap[operation(u)][x]), depvars, derivweights, (x2i(s, u, x), x), u, wind_ufunc, indexmap) for d in ( let orders = derivweights.orders[x] setdiff(orders[isodd.(orders)], skip) end - )] for x in params(u, s)], init = []) for u in depvars], init = []) + )] for x in ivs(u, s)], init = []) for u in depvars], init = []) ) wind_rules = [] @@ -106,7 +106,7 @@ end end return safe_vcat(wind_rules, vec(mapreduce(safe_vcat, depvars, init = []) do u - mapreduce(safe_vcat, params(u, s), init = []) do x + mapreduce(safe_vcat, ivs(u, s), init = []) do x j = x2i(s, u, x) let orders = setdiff(derivweights.orders[x], skip) oddorders = orders[isodd.(orders)] diff --git a/src/discretization/schemes/upwind_difference/upwind_difference_array.jl b/src/discretization/schemes/upwind_difference/upwind_difference_array.jl new file mode 100644 index 000000000..6e6906607 --- /dev/null +++ b/src/discretization/schemes/upwind_difference/upwind_difference_array.jl @@ -0,0 +1,135 @@ +######################################################################################## +# Stencil interface +######################################################################################## + +function _upwind_difference(D, ranges, interior, is, s, + bs, jx, u, udisc, ispositive) + args = ivs(u, s) + + j, x = jx + lenx = length(s, x) + haslower, hasupper = haslowerupper(bs, x) + + upperops = [] + lowerops = [] + if ispositive + if !haslower + lowerops = map(interior[j][1]:D.offside) do iboundary + lower_boundary_deriv(D, udisc, iboundary, j, is, interior) + end + end + interiorop = interior_deriv(D, udisc, s, -D.stencil_length+1:0, j, is, interior, bs) + else + if !hasupper + upperops = map((lenx-D.boundary_point_count+1):interior[j][end]) do iboundary + upper_boundary_deriv(D, udisc, iboundary, j, is, interior, lenx) + end + end + interiorop = interior_deriv(D, udisc, s, 0:D.stencil_length-1, j, is, interior, bs) + end + boundaryoppairs = safe_vcat(lowerops, upperops) + + NullBG_ArrayMaker(ranges, safe_vcat([Tuple(interior) => interiorop], boundaryoppairs))[interior...] +end + +""" +# upwind_difference +Generate a finite difference expression in `u` using the upwind difference at point `II::CartesianIndex` +in the direction of `x` +""" +function upwind_difference(d::Int, ranges, interior, is, s::DiscreteSpace, b, derivweights, + jx, u, udisc, ispositive) + j, x = jx + # return if this is an ODE + ndims(u, s) == 0 && return Fill(Num(0), ()) + D = if !ispositive + derivweights.windmap[1][Differential(x)^d] + else + derivweights.windmap[2][Differential(x)^d] + end + #@show D.stencil_coefs, D.stencil_length, D.boundary_stencil_length, D.boundary_point_count + # unit index in direction of the derivative + return _upwind_difference(D, ranges, interior, is, s, b, jx, u, udisc, ispositive) +end + +function upwind_difference(expr, d::Int, interior, s::DiscreteSpace, b, + depvars, derivweights, (j, x), u, udisc, indexmap) + # TODO: Allow derivatives in expr + + valrules = arrayvalmaps(s, u, depvars, interior) + exprarr = broadcast_substitute(expr, valrules) + ranges = get_ranges(u, s) + + IfElse.ifelse.(exprarr .> 0, + exprarr .* upwind_difference(d, ranges, get_interior(u, s, interior), get_is(u, s), s, b, derivweights, (j, x), u, udisc, true), + exprarr .* upwind_difference(d, ranges, get_interior(u, s, interior), get_is(u, s), s, b, derivweights, (j, x), u, udisc, false)) +end + +@inline function generate_winding_rules(interior, s::DiscreteSpace, depvars, + derivweights::DifferentialDiscretizer, bcmap, + indexmap, terms, skip = []) + # for all independent variables and dependant variables + rules = safe_vcat(#Catch multiplication + reduce(safe_vcat, + [reduce(safe_vcat, + [[@rule *(~~a, $(Differential(x)^d)(u), ~~b) => + upwind_difference(*(~a..., ~b...), d, interior, s, + bcmap[operation(u)][x], depvars, + derivweights, (x2i(s, u, x), x), u, + s.discvars[u], indexmap) + for d in (let orders = derivweights.orders[x] + setdiff(orders[isodd.(orders)], skip) + end)] + for x in ivs(u, s)], init = []) + for u in depvars], init = []), + + #Catch division and multiplication, see issue #1 + reduce(safe_vcat, + [reduce(safe_vcat, + [[@rule /(*(~~a, $(Differential(x)^d)(u), ~~b), ~c) => + upwind_difference(*(~a..., ~b...) / ~c, d, interior, s, + bcmap[operation(u)][x], depvars, + derivweights, (x2i(s, u, x), x), u, + s.discvars[u], indexmap) + for d in (let orders = derivweights.orders[x] + setdiff(orders[isodd.(orders)], skip) + + end)] + for x in ivs(u, s)], init = []) + for u in depvars], init = []) + ) + + wind_rules = [] + + # wind_exprs = [] + for t in terms + for r in rules + if r(t) !== nothing + push!(wind_rules, t => r(t)) + end + end + end + + return safe_vcat(wind_rules, vec(mapreduce(safe_vcat, depvars) do u + mapreduce(safe_vcat, ivs(u, s), init = []) do x + j = x2i(s, u, x) + is = get_is(u, s) + uinterior = get_interior(u, s, interior) + uranges = get_ranges(u, s) + let orders = derivweights.orders[x] + oddorders = setdiff(orders[isodd.(orders)], skip) + + # for all odd orders + if length(oddorders) > 0 + map(oddorders) do d + (Differential(x)^d)(u) => + upwind_difference(d, uranges, uinterior, is, s, bcmap[operation(u)][x], + derivweights, (j, x), u, s.discvars[u], true) + end + else + [] + end + end + end + end)) +end diff --git a/src/interface/scheme_types.jl b/src/interface/scheme_types.jl index 3d7c9035f..be8875ef0 100644 --- a/src/interface/scheme_types.jl +++ b/src/interface/scheme_types.jl @@ -27,3 +27,52 @@ function extent(::WENOScheme, dorder) end # Note: This type and its subtypes will become important later with the stencil interfaces as we will need to dispatch on derivative order and approximation order + +# Functional Schemes + +""" + # FunctionalScheme +A user definable scheme that takes a set of functions as input. The functions define the derivative at the interior, lower boundary, and upper boundary, taking the following inputs: + +Functions must be of the form `f(u, p, t, deriv_iv, d_iv)`. + +For the interior, `u` takes a vector of dependent variable values in the direction of the derivative +of length `interior_points`. `interior_points` must be odd, as this function defines the derivative at the center of the input points. + +For the lower and upper boundaries, `u` takes a vector of dependent variable values of length `boundary_points`. This will be the `boundary_points` number of points closest to the lower and upper boundary respectively. The boundary functions define the derivative at their index in the function vector, numbering from the boundary. For example, if `boundary_points = 3`, the first function in the vector will define the derivative at the boundary, the second at the boundary plus one step, and the third at the boundary plus two steps. +In general, `upper` and `lower` must be at least `floor(interior_points/2)` long. Where you have no good approximation for a derivative at the boundary, you can use `nothing` as a placeholder. MethodOfLines will then attempt to use an extrapolation here where nessesary. + +`p` will take all parameter values in the order specified in the PDESystem. + +`deriv_iv` takes a vector of independent variable values of the same support as for `u`, for the independent variable in the direction of the derivative. + +If `is_nonuniform` is false, `d_iv` will take a scalar value of the stepsize between the points used to call the function in `u` and `deriv_iv`. + +If `is_nonuniform` is true, the scheme must be able to accept `d_iv` as a vector of stepsizes between the points used to call the function in `u` and `deriv_iv`, therefore of length `length(u)-1`. If it is used on a uniform grid, the stepsizes will be the same for all points. +""" +struct FunctionalScheme{F, V1, V2} <: AbstractScheme + """ + The function that defines the scheme on the interior. + """ + interior::F + """ + The vector of functions that defines the scheme near the lower boundary. + """ + lower::V1 + """ + The vector functions that defines the scheme near the upper boundary. + """ + upper::V2 + """ + The number of interior points that the interior scheme takes as input + """ + interior_points::Int + """ + The number of boundary points that the lower and upper schemes take as input. + """ + boundary_points::Int + """ + Whether this scheme takes grid steps as input for the nonuniform case. + """ + is_nonuniform::Bool +end diff --git a/src/system_parsing/bcs/parse_boundaries.jl b/src/system_parsing/bcs/parse_boundaries.jl index db451dfb3..dbaac9b17 100644 --- a/src/system_parsing/bcs/parse_boundaries.jl +++ b/src/system_parsing/bcs/parse_boundaries.jl @@ -179,9 +179,9 @@ function generate_boundary_matching_rules(v, orders) upperboundary(x) = v.intervals[x][2] # Rules to match boundary conditions on the lower boundaries - lower = Dict([operation(u) => Dict([x => _boundary_rules(v, orders, u, x, lowerboundary(x)) for x in all_params(u, v)]) for u in v.ū]) + lower = Dict([operation(u) => Dict([x => _boundary_rules(v, orders, u, x, lowerboundary(x)) for x in all_ivs(u, v)]) for u in v.ū]) - upper = Dict([operation(u) => Dict([x => _boundary_rules(v, orders, u, x, upperboundary(x)) for x in all_params(u, v)]) for u in v.ū]) + upper = Dict([operation(u) => Dict([x => _boundary_rules(v, orders, u, x, upperboundary(x)) for x in all_ivs(u, v)]) for u in v.ū]) return (lower, upper) end diff --git a/src/system_parsing/variable_map.jl b/src/system_parsing/variable_map.jl index 7ec4a0646..4d8bdb573 100644 --- a/src/system_parsing/variable_map.jl +++ b/src/system_parsing/variable_map.jl @@ -1,6 +1,7 @@ struct VariableMap ū x̄ + ps time intervals args @@ -9,8 +10,11 @@ struct VariableMap i2x end -function VariableMap(eqs, depvars, domain, time) +function VariableMap(eqs, depvars, domain, time, ps = []) time = safe_unwrap(time) + ps = map(ps) do p + safe_unwrap(p.first) + end depvar_ops = get_ops(depvars) # Get all dependent variables in the correct type alldepvars = get_all_depvars(eqs, depvar_ops) @@ -27,7 +31,7 @@ function VariableMap(eqs, depvars, domain, time) args = [operation(u) => arguments(u) for u in ū] x̄2dim = [x̄[i] => i for i in 1:nspace] dim2x̄ = [i => x̄[i] for i in 1:nspace] - return VariableMap(ū, x̄, time, Dict(intervals), Dict(args), depvar_ops, Dict(x̄2dim), Dict(dim2x̄)) + return VariableMap(ū, x̄, ps, time, Dict(intervals), Dict(args), depvar_ops, Dict(x̄2dim), Dict(dim2x̄)) end VariableMap(pdesys::PDESystem, disc::MOLFiniteDifference) = VariableMap(pdesys.eqs, pdesys.dvs, pdesys.domain, disc.time) @@ -40,13 +44,13 @@ function update_varmap!(v, newdv) end -params(u, v::VariableMap) = remove(v.args[operation(u)], v.time) +ivs(u, v::VariableMap) = remove(v.args[operation(u)], v.time) -Base.ndims(u, v::VariableMap) = length(params(u, v)) +Base.ndims(u, v::VariableMap) = length(ivs(u, v)) all_ivs(v::VariableMap) = v.time === nothing ? v.x̄ : v.x̄ ∪ [v.time] -all_params(u, v::VariableMap) = v.args[operation(u)] +all_ivs(u, v::VariableMap) = v.args[operation(u)] depvar(u, v::VariableMap) = operation(u)(v.args[operation(u)]...) @@ -54,7 +58,7 @@ x2i(v::VariableMap, u, x) = findfirst(isequal(x), remove(v.args[operation(u)], v @inline function axiesvals(v::VariableMap, u_, x_, I) u = depvar(u_, v) - map(params(u, v)) do x + map(ivs(u, v)) do x x => (I[x2i(v, u, x)] == 1 ? v.intervals[x][1] : v.intervals[x][2]) end end diff --git a/test/components/finite_diff_schemes.jl b/test/components/finite_diff_schemes.jl index 256d62b14..66fb5f5e7 100644 --- a/test/components/finite_diff_schemes.jl +++ b/test/components/finite_diff_schemes.jl @@ -1,6 +1,6 @@ using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils, LinearAlgebra -@parameters x, t +@parameters x, t @variables u(..) Dx(d) = Differential(x)^d @@ -10,7 +10,7 @@ Dt = Differential(t) t_min= 0. t_max = 2.0 x_min = 0. -x_max = 20.0 +x_max = 20.0 dx = 1.0 @@ -28,7 +28,7 @@ domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)] @named pdesys = PDESystem(pde,bcs,domains,[t,x],[u(t,x)]) - # Test centered order + # Test centered order disc = MOLFiniteDifference([x=>dx], t; approx_order=a) depvar_ops = map(x->operation(x.val),pdesys.depvars) @@ -44,7 +44,7 @@ domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)] s = MethodOfLines.DiscreteSpace(domains, depvars, x̄, disc) derivweights = MethodOfLines.DifferentialDiscretizer(pdesys, s, disc) - + #@show pde.rhs, operation(pde.rhs), arguments(pde.rhs) for II in s.Igrid[s.ū[1]][2:end-1] #II = s.Igrid[end-1] @@ -67,16 +67,16 @@ domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)] end @testset "Test 01: Nonlinear Diffusion discretization" begin - + pde = Dt(u(t,x)) ~ Dx(1)(u(t,x)) bcs = [u(0,x) ~ cos(x), u(t,0) ~ exp(-t), u(t,Float64(π)) ~ -exp(-t)] @named pdesys = PDESystem(pde,bcs,domains,[t,x],[u(t,x)]) - # Test centered order - + # Test centered order + depvar_ops = map(x->operation(x.val),pdesys.depvars) - + depvars_lhs = MethodOfLines.get_depvars(pde.lhs, depvar_ops) depvars_rhs = MethodOfLines.get_depvars(pde.rhs, depvar_ops) depvars = collect(depvars_lhs ∪ depvars_rhs) @@ -84,15 +84,15 @@ end # ignore if the only argument is [t] indvars = first(Set(filter(xs->!isequal(xs, [t]), map(arguments, depvars)))) x̄ = first(Set(filter(!isempty, map(u->filter(x-> t === nothing || !isequal(x, t.val), arguments(u)), depvars)))) - + for order in [2] disc = MOLFiniteDifference([x=>dx], t; approx_order=order) s = MethodOfLines.DiscreteSpace(domains, depvars, x̄, disc) derivweights = MethodOfLines.DifferentialDiscretizer(pdesys, s, disc) - + ufunc(u, I, x) = s.discvars[u][I] - #TODO Test Interpolation of params + #TODO Test Interpolation of ivs # Test simple case for II in s.Igrid[s.ū[1]][2:end-1] expr = MethodOfLines.cartesian_nonlinear_laplacian((1~1).lhs, II, derivweights, s, depvars, x, u(t,x)) @@ -103,14 +103,14 @@ end end @testset "Test 02: Spherical Diffusion discretization" begin - + pde = Dt(u(t,x)) ~ 1/x^2 * Dx(1)(x^2 * Dx(1)(u(t,x))) bcs = [u(0,x) ~ cos(x), u(t,0) ~ exp(-t), u(t,Float64(π)) ~ -exp(-t)] @named pdesys = PDESystem(pde,bcs,domains,[t,x],[u(t,x)]) - # Test centered order + # Test centered order disc = MOLFiniteDifference([x=>dx], t; approx_order=2) depvar_ops = map(x->operation(x.val),pdesys.depvars) @@ -126,11 +126,11 @@ end s = MethodOfLines.DiscreteSpace(domains, depvars, x̄, disc) derivweights = MethodOfLines.DifferentialDiscretizer(pdesys, s, disc) - + for II in s.Igrid[s.ū[1]][2:end-1] - #TODO Test Interpolation of params + #TODO Test Interpolation of ivs expr = MethodOfLines.spherical_diffusion((1~1).lhs, II, derivweights, s, depvars, x, u(t,x)) #@show II, expr end -end \ No newline at end of file +end diff --git a/test/components/function_scheme.jl b/test/components/function_scheme.jl new file mode 100644 index 000000000..fc5c1facb --- /dev/null +++ b/test/components/function_scheme.jl @@ -0,0 +1,11 @@ +@testset "User defined scheme, 3 interior points" begin + f_interior = (u, p, t, x, dx) -> IfElse.ifelse(u[2] < 0, (u[2] - u[1]) / dx, (u[3] - u[2]) / dx) + f_lower = (u, p, t, x, dx) -> (u[2] - u[1]) / dx + f_upper = (u, p, t, x, dx) -> (u[2] - u[1]) / dx + + scheme = FunctionScheme{3, 2}(f_interior, f_lower, f_upper, false) + + disc = MOLFiniteDifference([x => dx], t, advection_scheme=scheme) + + +end From 93fd03ac5aa10b22857986de78ad1739d450584f Mon Sep 17 00:00:00 2001 From: xtalax Date: Thu, 2 Mar 2023 20:25:58 +0000 Subject: [PATCH 02/26] remove dead files --- .../array_form/stencil_utils.jl | 174 ------------------ .../generate_finite_difference_rules.jl | 8 +- src/discretization/schemes/WENO/WENO.jl | 55 ++---- src/discretization/schemes/WENO/WENO_array.jl | 89 --------- .../function_scheme/function_scheme.jl | 56 ++++-- .../integral_expansion_array.jl | 46 ----- .../upwind_difference_array.jl | 135 -------------- src/interface/scheme_types.jl | 74 +++++--- 8 files changed, 105 insertions(+), 532 deletions(-) delete mode 100644 src/discretization/array_form/stencil_utils.jl delete mode 100644 src/discretization/schemes/WENO/WENO_array.jl delete mode 100644 src/discretization/schemes/integral_expansion/integral_expansion_array.jl delete mode 100644 src/discretization/schemes/upwind_difference/upwind_difference_array.jl diff --git a/src/discretization/array_form/stencil_utils.jl b/src/discretization/array_form/stencil_utils.jl deleted file mode 100644 index 7009dcce2..000000000 --- a/src/discretization/array_form/stencil_utils.jl +++ /dev/null @@ -1,174 +0,0 @@ -function lower_boundary_deriv(D, udisc, iboundary, j, is, interior) - weights = D.low_boundary_coefs[iboundary] - taps = 1:D.boundary_stencil_length - prepare_boundary_op((BoundaryDerivArrayOp(weights, taps, udisc, j, is, interior), - iboundary), interior, j) -end - -function upper_boundary_deriv(D, udisc, iboundary, j, is, interior, lenx) - weights = D.high_boundary_coefs[lenx - iboundary + 1] - taps = (lenx-D.boundary_stencil_length+1):lenx - prepare_boundary_op((BoundaryDerivArrayOp(weights, taps, udisc, j, is, interior), - iboundary), interior, j) -end - -function integral_op_pair(dx, udisc, j, is, ranges, interior, i) - prepare_boundary_op((IntegralArrayOp(dx, udisc, i, j, is, ranges, interior), - i), interior, j) -end - -function prepare_boundary_op(boundaryop, interior, j) - function maketuple(i) - out = map(1:length(interior)) do k - k == j ? i : interior[k] - end - return Tuple(out) - end - (op, iboundary) = boundaryop - return maketuple(iboundary) => op -end - -function prepare_boundary_ops(boundaryops, interior, j) - function maketuple(i) - out = map(1:length(interior)) do k - k == j ? i : interior[k] - end - return Tuple(out) - end - return map(boundaryops) do (op, iboundary) - maketuple(iboundary) => op - end -end - -function interior_deriv(D::DerivativeOperator{T,N,Wind,DX}, udisc, s, offsets, j, is, interior, bs, isx = false) where {T,N,Wind,DX<:Number} - weights = D.stencil_coefs - taps = offsets .+ is[j] - InteriorDerivArrayOp(weights, taps, udisc, s, j, is, interior, bs, isx) -end - -function interior_deriv(D::DerivativeOperator{T,N,Wind,DX}, udisc, s, offsets, j, is, interior, bs, isx = false) where {T,N,Wind,DX<:AbstractVector} - @assert length(bs) == 0 "Interface boundary conditions are not yet supported for nonuniform dx dimensions, such as $x, please post an issue to https://github.com/SciML/MethodOfLines.jl if you need this functionality." - weights = D.stencil_coefs[is[j]-D.boundary_point_count] - taps = offsets .+ is[j] - InteriorDerivArrayOp(weights, taps, udisc, s, j, is, interior, bs, isx) -end - -function BoundaryDerivArrayOp(weights, taps, udisc, j, is, interior) - # * I Possibly needs updating - Is = map(taps) do tap - map(1:ndims(udisc)) do i - if i == j - tap - else - is[i] - end - end - end - expr = sym_dot(weights, map(I -> udisc[I...], Is)) - - symindices = setdiff(1:ndims(udisc), [j]) - output_idx = Tuple(is[symindices]) - return FillArrayOp(recursive_unwrap(expr), output_idx, interior[symindices]) -end - -reduce_interior(interior, j) = first(interior[j]) == 1 ? [interior[1:j-1]..., 2:last(interior[j]), interior[j+1:end]...] : interior - -function trapezium_sum(ranges, interior, dx, udisc, is, j) - N = ndims(udisc) - I1 = unitindex(N, j) - I = CartesianIndex(is...) - i = is[j] - Im1 = I - I1 - im1 = is[j] - 1 - rinterior = reduce_interior(interior, j) - if dx isa Number - expr = (dx*(udisc[Im1] + udisc[I]) / 2) - else - expr = (dx[im1]*udisc[Im1] + dx[i]*udisc[I]) / 2 - end - - return FillArrayMaker(expr, is, ranges, rinterior)[interior...] -end - -function IntegralArrayOp(dx, udisc, k, j, is, ranges, interior, iswd = false) - if iswd - interior = [interior[1:j-1]..., 1:size(udisc, j), interior[j+1:end]...] - end - trapop = trapezium_sum(ranges, interior, dx, udisc, is, j) - op = sum(trapop[1:(k-first(interior[j])+1)], dims = j) - - return op -end - -function IntegralArrayMaker(dx, udisc, k, j, is, ranges, interior, iswd = false, ) - if iswd - ranges = [ranges[1:j-1]..., 1:size(udisc, j), ranges[j+1:end]...] - end - - op = IntegralArrayOp(dx, udisc, k, j, is, ranges, interior, iswd) - return FillArrayMaker(op, is, ranges, interior) -end - -function InteriorDerivArrayOp(weights, taps, udisc, s, j, output_idx, interior, bs, isx = false) - # * I Possibly needs updating - Is = map(taps) do tap - _is = map(1:ndims(udisc)) do i - if i == j - tap - else - output_idx[i] - end - end - CartesianIndex(_is...) - end - # Wrap interfaces - Is = map(I -> bwrap(I, bs, s, j, isx), Is) - - Is = map(Is) do I - map(1:ndims(udisc)) do i - I[i] - end - end - - expr = sym_dot(weights, map(I -> udisc[I...], Is)) - - return FillArrayOp(expr, output_idx, interior) -end - -function FillArrayOp(expr, output_idx, interior) - ranges = Dict(output_idx .=> interior) # hope this doesn't check bounds eagerly - return ArrayOp(Array{symtype(expr),length(output_idx)}, - output_idx, expr, +, nothing, ranges) -end - -NullBG_ArrayMaker(ranges, ops) = ArrayMaker{Real}(Tuple(map(r -> r[end] - r[1] + 1, ranges)), - vcat(Tuple(ranges) => 0, ops)) -Construct_ArrayMaker(ranges, - ops) = ArrayMaker{Real}(Tuple(map(r -> r[end] - r[1] + 1, ranges)), ops) -Construct_ArrayMaker{T}(ranges, - ops) = ArrayMaker{T}(Tuple(map(r -> r[end] - r[1] + 1, ranges)), ops) - - - FillArrayMaker(expr, is, - ranges, interior) = NullBG_ArrayMaker(ranges, - [Tuple(interior) => FillArrayOp(expr, is, interior)]) - -ArrayMakerWrap(udisc, ranges) = Arraymaker{Real}(Tuple(map(r -> r[end] - r[1] + 1, ranges)), - [Tuple(ranges) => udisc]) - -##### - -function get_interior(u, s, interior) - map(ivs(u, s)) do x - if haskey(interior, x) - interior[x] - else - 1:length(s, x) - end - end -end - -function get_ranges(u, s) - map(x -> first(axes(s.grid[x])), ivs(u, s)) -end -get_is(u, s) = map(x -> s.index_syms[x], ivs(u, s)) diff --git a/src/discretization/generate_finite_difference_rules.jl b/src/discretization/generate_finite_difference_rules.jl index 379be5a89..e443b83b3 100644 --- a/src/discretization/generate_finite_difference_rules.jl +++ b/src/discretization/generate_finite_difference_rules.jl @@ -54,10 +54,12 @@ function generate_finite_difference_rules(II::CartesianIndex, s::DiscreteSpace, # Advection rules if derivweights.advection_scheme isa UpwindScheme advection_rules = generate_winding_rules(II, s, depvars, derivweights, bmap, indexmap, terms) - elseif derivweights.advection_scheme isa WENOScheme - advection_rules = generate_WENO_rules(II, s, depvars, derivweights, bmap, indexmap, terms) + elseif derivweights.advection_scheme isa FunctionalScheme + advection_rules = generate_advection_rules(derivweights.advection_scheme, II, s, + depvars, derivweights, bmap, indexmap, terms) advection_rules = vcat(advection_rules, - generate_winding_rules(II, s, depvars, derivweights, bmap, indexmap, terms; skip = [1])) + generate_winding_rules(II, s, depvars, derivweights, bmap, + indexmap, terms; skip = [1])) else error("Unsupported advection scheme $(derivweights.advection_scheme) encountered.") end diff --git a/src/discretization/schemes/WENO/WENO.jl b/src/discretization/schemes/WENO/WENO.jl index 827d0e72f..ff62f383b 100644 --- a/src/discretization/schemes/WENO/WENO.jl +++ b/src/discretization/schemes/WENO/WENO.jl @@ -1,37 +1,11 @@ +function weno_f(u, p, t, x, dx) + ε = p[1] -""" -Implements the WENO scheme of Jiang and Shu. - -Specified in https://repository.library.brown.edu/studio/item/bdr:297524/PDF/ (Page 8-9) - -Implementation *heavily* inspired by https://github.com/ranocha/HyperbolicDiffEq.jl/blob/84c2d882e0c8956457c7d662bf7f18e3c27cfa3d/src/finite_volumes/weno_jiang_shu.jl by H. Ranocha. -""" -function weno(II::CartesianIndex, s::DiscreteSpace, wenoscheme::WENOScheme, bs, jx, u, dx::Number) - j, x = jx - ε = wenoscheme.epsilon - - I1 = unitindex(ndims(u, s), j) - - udisc = s.discvars[u] - - Im2 = bwrap(II - 2I1, bs, s, jx) - Im1 = bwrap(II - I1, bs, s, jx) - Ip1 = bwrap(II + I1, bs, s, jx) - Ip2 = bwrap(II + 2I1, bs, s, jx) - is = map(I -> I[j], [Im2, Im1, Ip1, Ip2]) - for i in is - if i < 1 - return nothing - elseif i > length(s, x) - return nothing - end - end - - u_m2 = udisc[Im2] - u_m1 = udisc[Im1] - u_0 = udisc[II] - u_p1 = udisc[Ip1] - u_p2 = udisc[Ip2] + u_m2 = u[1] + u_m1 = u[2] + u_0 = u[3] + u_p1 = u[4] + u_p2 = u[5] γm1 = 1 / 10 γm2 = 3 / 5 @@ -74,16 +48,15 @@ function weno(II::CartesianIndex, s::DiscreteSpace, wenoscheme::WENOScheme, bs, hp = wp1 * hp1 + wp2 * hp2 + wp3 * hp3 hm = wm1 * hm1 + wm2 * hm2 + wm3 * hm3 - return recursive_unwrap((hp - hm) / dx) -end - -function weno(II::CartesianIndex, s::DiscreteSpace, b, jx, u, dx::AbstractVector) - @assert false "WENO scheme not implemented for nonuniform grids." + return (hp - hm) / dx end """ -This is a catch all ruleset, as such it does not use @rule. +`WENOScheme` of Jiang and Shu +## Keyword Arguments +- `epsilon`: A quantity used to prevent vanishing denominators in the scheme, defaults to `1e-6`. More sensetive problems will benefit from a smaller value. It is defined as a functional scheme. """ -@inline function generate_WENO_rules(II::CartesianIndex, s::DiscreteSpace, depvars, derivweights::DifferentialDiscretizer, bcmap, indexmap, terms) - return reduce(safe_vcat, [[(Differential(x))(u) => weno(Idx(II, s, u, indexmap), s, derivweights.advection_scheme, filter_interfaces(bcmap[operation(u)][x]), (x2i(s, u, x), x), u, s.dxs[x]) for x in ivs(u, s)] for u in depvars], init = []) +function WENOScheme(epsilon = 1e-6) + boundry_f = [nothing, nothing] + return FunctionScheme{5, 0}(weno_f, boundary_f, boundary_f, false, [epsilon], name = "WENO") end diff --git a/src/discretization/schemes/WENO/WENO_array.jl b/src/discretization/schemes/WENO/WENO_array.jl deleted file mode 100644 index 31012fcc4..000000000 --- a/src/discretization/schemes/WENO/WENO_array.jl +++ /dev/null @@ -1,89 +0,0 @@ -######################################################################################## -# Stencil interface -######################################################################################## - -function weno(interior, s::DiscreteSpace, wenoscheme::WENOScheme, bs, jx, u, dx::Number) - j, x = jx - ε = wenoscheme.epsilon - - interior = get_interior(u, s, interior) - ranges = get_ranges(u, s) - is = get_is(u, s) - - II = CartesianIndex(is...) - I1 = unitindex(ndims(u, s), j) - - udisc = s.discvars[u] - - Im2 = bwrap(II - 2I1, bs, s, j) - Im1 = bwrap(II - I1, bs, s, j) - Ip1 = bwrap(II + I1, bs, s, j) - Ip2 = bwrap(II + 2I1, bs, s, j) - - u_m2 = udisc[Im2] - u_m1 = udisc[Im1] - u_0 = udisc[II] - u_p1 = udisc[Ip1] - u_p2 = udisc[Ip2] - - γm1 = 1 / 10 - γm2 = 3 / 5 - γm3 = 3 / 10 - - β1 = 13 * (u_0 - 2 * u_p1 + u_p2)^2 / 12 + (3 * u_0 - 4 * u_p1 + u_p2)^2 / 4 - β2 = 13 * (u_m1 - 2 * u_0 + u_p1)^2 / 12 + (u_m1 - u_p1)^2 / 4 - β3 = 13 * (u_m2 - 2 * u_m1 + u_0)^2 / 12 + (u_m2 - 4 * u_m1 + 3 * u_0)^2 / 4 - - ωm1 = γm1 / (ε + β1)^2 - ωm2 = γm2 / (ε + β2)^2 - ωm3 = γm3 / (ε + β3)^2 - - wm_denom = ωm1 + ωm2 + ωm3 - wm1 = ωm1 / wm_denom - wm2 = ωm2 / wm_denom - wm3 = ωm3 / wm_denom - - γp1 = 3 / 10 - γp2 = 3 / 5 - γp3 = 1 / 10 - - ωp1 = γp1 / (ε + β1)^2 - ωp2 = γp2 / (ε + β2)^2 - ωp3 = γp3 / (ε + β3)^2 - - wp_denom = ωp1 + ωp2 + ωp3 - wp1 = ωp1 / wp_denom - wp2 = ωp2 / wp_denom - wp3 = ωp3 / wp_denom - - hm1 = (11 * u_0 - 7 * u_p1 + 2 * u_p2) / 6 - hm2 = (5 * u_0 - u_p1 + 2 * u_m1) / 6 - hm3 = (2 * u_0 + 5 * u_m1 - u_m2) / 6 - - hp1 = (2 * u_0 + 5 * u_p1 - u_p2) / 6 - hp2 = (5 * u_0 + 2 * u_p1 - u_m1) / 6 - hp3 = (11 * u_0 - 7 * u_m1 + 2 * u_m2) / 6 - - hp = wp1 * hp1 + wp2 * hp2 + wp3 * hp3 - hm = wm1 * hm1 + wm2 * hm2 + wm3 * hm3 - - expr = (hp - hm) / dx - return FillArrayMaker(recursive_unwrap(expr), is, ranges, interior)[interior...] -end - -function weno(interior, s::DiscreteSpace, b, jx, u, dx::AbstractVector) - @assert false "WENO scheme not implemented for nonuniform grids." -end - -""" -This is a catch all ruleset, as such it does not use @rule. -""" -@inline function generate_WENO_rules(interior, s::DiscreteSpace, depvars, - derivweights::DifferentialDiscretizer, pmap, - indexmap, terms) - return reduce(safe_vcat, - [[(Differential(x))(u) => weno(interior, s, derivweights.advection_scheme, - pmap.map[operation(u)][x], (x2i(s, u, x), x), u, s.dxs[x]) - for x in ivs(u, s)] - for u in depvars], init = []) -end diff --git a/src/discretization/schemes/function_scheme/function_scheme.jl b/src/discretization/schemes/function_scheme/function_scheme.jl index b56078c1f..d048ac9b9 100644 --- a/src/discretization/schemes/function_scheme/function_scheme.jl +++ b/src/discretization/schemes/function_scheme/function_scheme.jl @@ -1,19 +1,5 @@ -function function_scheme(F::FunctionScheme{is_nonuniform}, II, s, bs, jx, u, ufunc) where {is_nonuniform <: Val{false}} - - - # Tap points of the stencil, this uses boundary_point_count as this is equal to half the stencil size, which is what we want. - u_disc = ufunc(u, Itap, x) - ps = params(s) - dx = s.dxs[x] - t = s.time - - return f(u_disc, ps, t, x, dx) -end - - -function get_f_and_taps(F::FunctionScheme{is_nonuniform}, II, s, bs, jx, u, ufunc) where {is_nonuniform <: Val{false}} +function get_f_and_taps(F::FunctionalScheme, II, s, bs, jx, u) j, x = jx - ndims(u, s) == 0 && return 0 # unit index in direction of the derivative I1 = unitindex(ndims(u, s), j) # offset is important due to boundary proximity @@ -22,7 +8,7 @@ function get_f_and_taps(F::FunctionScheme{is_nonuniform}, II, s, bs, jx, u, ufun lower_point_count = length(f.lower) upper_point_count = length(f.upper) - if (II[j] <= D.lower_point_count) & !haslower + if (II[j] <= lower_point_count) & !haslower f = F.lower[II[j]] offset = 1 - II[j] Itap = [II + (i + offset) * I1 for i in 0:(F.boundary_points-1)] @@ -36,3 +22,41 @@ function get_f_and_taps(F::FunctionScheme{is_nonuniform}, II, s, bs, jx, u, ufun end return f, Itap +end + +function function_scheme(F::FunctionalScheme{is_nonuniform}, II, s, bs, jx, u, ufunc) where {is_nonuniform<:Val{false}} + j, x = jx + ndims(u, s) == 0 && return 0 + + f, Itap = get_f_and_taps(F, II, s, bs, jx, u) + if isnothing(f) + error("Scheme $(F.name) applied to $u in direction of $x at point $II is not defined.") + end + # Tap points of the stencil, this uses boundary_point_count as this is equal to half the stencil size, which is what we want. + u_disc = ufunc(u, Itap, x) + ps = params(s) + dx = s.dxs[x] + if F.is_nonuniform + if dx isa AbstractVector + dx = map(I -> dx[I[j]], Itap)# + end + elseif dx isa AbstractVector + error("Scheme $(F.name) not implemented for nonuniform dxs.") + end + + t = s.time + + return f(u_disc, ps, t, x, dx) +end + +@inline function generate_advection_rules(F::FunctionalScheme, II::CartesianIndex, s::DiscreteSpace, depvars, derivweights::DifferentialDiscretizer, bcmap, indexmap, terms) + central_ufunc(u, I, x) = s.discvars[u][I] + return reduce(safe_vcat, [[(Differential(x))(u) => + function_scheme(F, + Idx(II, s, u, indexmap), s, + filter_interfaces(bcmap[operation(u)][x]), + (x2i(s, u, x), x), u, + central_ufunc) + for x in ivs(u, s)] + for u in depvars], init=[]) +end diff --git a/src/discretization/schemes/integral_expansion/integral_expansion_array.jl b/src/discretization/schemes/integral_expansion/integral_expansion_array.jl deleted file mode 100644 index d10f5b034..000000000 --- a/src/discretization/schemes/integral_expansion/integral_expansion_array.jl +++ /dev/null @@ -1,46 +0,0 @@ -function euler_integral(interior, s, jx, u, udisc) - j, x = jx - dx = s.dxs[x] - interior = get_interior(u, s, interior) - ranges = get_ranges(u, s) - is = get_is(u, s) - - oppairs = map(interior[j]) do i - integral_op_pair(dx, udisc, j, is, ranges, interior, i) - end - - return NullBG_ArrayMaker(ranges, oppairs)[interior...] -end - -# An integral across the whole domain (xmin .. xmax) -function whole_domain_integral(interior, s, jx, u, udisc) - j, x = jx - dx = s.dxs[x] - - interior = get_interior(u, s, interior) - ranges = get_ranges(u, s) - is = get_is(u, s) - - lenx = length(s, x) - - return IntegralArrayMaker(dx, udisc, lenx, j, is, ranges, interior, true)[interior...] -end - -@inline function generate_euler_integration_rules(interior, s::DiscreteSpace, depvars, indexmap, terms) - eulerrules = reduce(safe_vcat, [[Integral(x in DomainSets.ClosedInterval(s.vars.intervals[x][1], Num(x)))(u) => - euler_integral(interior, s, (x2i(s, u, x), x), u, s.discvars[u]) - for x in ivs(u, s)] - for u in depvars], init = []) - return eulerrules -end - -@inline function generate_whole_domain_integration_rules(interior, s::DiscreteSpace, depvars, indexmap, terms, bvar = nothing) - wholedomainrules = reduce(safe_vcat, - [[Integral(x in DomainSets.ClosedInterval(s.vars.intervals[x][1], s.vars.intervals[x][2]))(u) => - whole_domain_integral(interior, s, (x2i(s, u, x), x), u, s.discvars[u]) - for x in filter(x -> (!haskey(indexmap, x) | isequal(x, bvar)), ivs(u, s))] - for u in depvars], - init = []) - - return wholedomainrules -end diff --git a/src/discretization/schemes/upwind_difference/upwind_difference_array.jl b/src/discretization/schemes/upwind_difference/upwind_difference_array.jl deleted file mode 100644 index 6e6906607..000000000 --- a/src/discretization/schemes/upwind_difference/upwind_difference_array.jl +++ /dev/null @@ -1,135 +0,0 @@ -######################################################################################## -# Stencil interface -######################################################################################## - -function _upwind_difference(D, ranges, interior, is, s, - bs, jx, u, udisc, ispositive) - args = ivs(u, s) - - j, x = jx - lenx = length(s, x) - haslower, hasupper = haslowerupper(bs, x) - - upperops = [] - lowerops = [] - if ispositive - if !haslower - lowerops = map(interior[j][1]:D.offside) do iboundary - lower_boundary_deriv(D, udisc, iboundary, j, is, interior) - end - end - interiorop = interior_deriv(D, udisc, s, -D.stencil_length+1:0, j, is, interior, bs) - else - if !hasupper - upperops = map((lenx-D.boundary_point_count+1):interior[j][end]) do iboundary - upper_boundary_deriv(D, udisc, iboundary, j, is, interior, lenx) - end - end - interiorop = interior_deriv(D, udisc, s, 0:D.stencil_length-1, j, is, interior, bs) - end - boundaryoppairs = safe_vcat(lowerops, upperops) - - NullBG_ArrayMaker(ranges, safe_vcat([Tuple(interior) => interiorop], boundaryoppairs))[interior...] -end - -""" -# upwind_difference -Generate a finite difference expression in `u` using the upwind difference at point `II::CartesianIndex` -in the direction of `x` -""" -function upwind_difference(d::Int, ranges, interior, is, s::DiscreteSpace, b, derivweights, - jx, u, udisc, ispositive) - j, x = jx - # return if this is an ODE - ndims(u, s) == 0 && return Fill(Num(0), ()) - D = if !ispositive - derivweights.windmap[1][Differential(x)^d] - else - derivweights.windmap[2][Differential(x)^d] - end - #@show D.stencil_coefs, D.stencil_length, D.boundary_stencil_length, D.boundary_point_count - # unit index in direction of the derivative - return _upwind_difference(D, ranges, interior, is, s, b, jx, u, udisc, ispositive) -end - -function upwind_difference(expr, d::Int, interior, s::DiscreteSpace, b, - depvars, derivweights, (j, x), u, udisc, indexmap) - # TODO: Allow derivatives in expr - - valrules = arrayvalmaps(s, u, depvars, interior) - exprarr = broadcast_substitute(expr, valrules) - ranges = get_ranges(u, s) - - IfElse.ifelse.(exprarr .> 0, - exprarr .* upwind_difference(d, ranges, get_interior(u, s, interior), get_is(u, s), s, b, derivweights, (j, x), u, udisc, true), - exprarr .* upwind_difference(d, ranges, get_interior(u, s, interior), get_is(u, s), s, b, derivweights, (j, x), u, udisc, false)) -end - -@inline function generate_winding_rules(interior, s::DiscreteSpace, depvars, - derivweights::DifferentialDiscretizer, bcmap, - indexmap, terms, skip = []) - # for all independent variables and dependant variables - rules = safe_vcat(#Catch multiplication - reduce(safe_vcat, - [reduce(safe_vcat, - [[@rule *(~~a, $(Differential(x)^d)(u), ~~b) => - upwind_difference(*(~a..., ~b...), d, interior, s, - bcmap[operation(u)][x], depvars, - derivweights, (x2i(s, u, x), x), u, - s.discvars[u], indexmap) - for d in (let orders = derivweights.orders[x] - setdiff(orders[isodd.(orders)], skip) - end)] - for x in ivs(u, s)], init = []) - for u in depvars], init = []), - - #Catch division and multiplication, see issue #1 - reduce(safe_vcat, - [reduce(safe_vcat, - [[@rule /(*(~~a, $(Differential(x)^d)(u), ~~b), ~c) => - upwind_difference(*(~a..., ~b...) / ~c, d, interior, s, - bcmap[operation(u)][x], depvars, - derivweights, (x2i(s, u, x), x), u, - s.discvars[u], indexmap) - for d in (let orders = derivweights.orders[x] - setdiff(orders[isodd.(orders)], skip) - - end)] - for x in ivs(u, s)], init = []) - for u in depvars], init = []) - ) - - wind_rules = [] - - # wind_exprs = [] - for t in terms - for r in rules - if r(t) !== nothing - push!(wind_rules, t => r(t)) - end - end - end - - return safe_vcat(wind_rules, vec(mapreduce(safe_vcat, depvars) do u - mapreduce(safe_vcat, ivs(u, s), init = []) do x - j = x2i(s, u, x) - is = get_is(u, s) - uinterior = get_interior(u, s, interior) - uranges = get_ranges(u, s) - let orders = derivweights.orders[x] - oddorders = setdiff(orders[isodd.(orders)], skip) - - # for all odd orders - if length(oddorders) > 0 - map(oddorders) do d - (Differential(x)^d)(u) => - upwind_difference(d, uranges, uinterior, is, s, bcmap[operation(u)][x], - derivweights, (j, x), u, s.discvars[u], true) - end - else - [] - end - end - end - end)) -end diff --git a/src/interface/scheme_types.jl b/src/interface/scheme_types.jl index be8875ef0..2a68e7e23 100644 --- a/src/interface/scheme_types.jl +++ b/src/interface/scheme_types.jl @@ -2,55 +2,40 @@ abstract type AbstractScheme end struct UpwindScheme <: AbstractScheme order - function UpwindScheme(approx_order = 1) + function UpwindScheme(approx_order=1) return new(approx_order) end end -extent(scheme::UpwindScheme, dorder) = dorder+scheme.order-1 - -""" -`WENOScheme` of Jiang and Shu -## Keyword Arguments -- `epsilon`: A quantity used to prevent vanishing denominators in the scheme, defaults to `1e-6`. More sensetive problems will benefit from a smaller value. -""" -struct WENOScheme <: AbstractScheme - epsilon - function WENOScheme(epsilon = 1e-6) - new(epsilon) - end -end - -function extent(::WENOScheme, dorder) - @assert dorder == 1 "Only first order spatial derivatives are implemented for WENO schemes." - return 2 -end - -# Note: This type and its subtypes will become important later with the stencil interfaces as we will need to dispatch on derivative order and approximation order +extent(scheme::UpwindScheme, dorder) = dorder + scheme.order - 1 # Functional Schemes """ # FunctionalScheme -A user definable scheme that takes a set of functions as input. The functions define the derivative at the interior, lower boundary, and upper boundary, taking the following inputs: +A user definable scheme that takes a set of functions as input. The functions define the derivative at the interior, lower boundary, and upper boundary. + +`lower` and `upper` should be vectors of functions. In general, `upper` and `lower` must be at least `floor(interior_points/2)` long. Where you have no good approximation for a derivative at the boundary, you can use `nothing` as a placeholder. MethodOfLines will then attempt to use an extrapolation here where nessesary. Be warned that this can lead to instability. + +The boundary functions define the derivative at their index in the function vector, numbering from the boundary. For example, if `boundary_points = 3`, the first function in the vector will define the derivative at the boundary, the second at the boundary plus one step, and the third at the boundary plus two steps. + +The functions making up the taking the following inputs: Functions must be of the form `f(u, p, t, deriv_iv, d_iv)`. For the interior, `u` takes a vector of dependent variable values in the direction of the derivative of length `interior_points`. `interior_points` must be odd, as this function defines the derivative at the center of the input points. -For the lower and upper boundaries, `u` takes a vector of dependent variable values of length `boundary_points`. This will be the `boundary_points` number of points closest to the lower and upper boundary respectively. The boundary functions define the derivative at their index in the function vector, numbering from the boundary. For example, if `boundary_points = 3`, the first function in the vector will define the derivative at the boundary, the second at the boundary plus one step, and the third at the boundary plus two steps. -In general, `upper` and `lower` must be at least `floor(interior_points/2)` long. Where you have no good approximation for a derivative at the boundary, you can use `nothing` as a placeholder. MethodOfLines will then attempt to use an extrapolation here where nessesary. - -`p` will take all parameter values in the order specified in the PDESystem. +For the lower and upper boundaries, `u` takes a vector of dependent variable values of length `boundary_points`. This will be the `boundary_points` number of points closest to the lower and upper boundary respectively. +`p` will take all parameter values in the order specified in the PDESystem, with the scheme's parameters prepended to the list. `deriv_iv` takes a vector of independent variable values of the same support as for `u`, for the independent variable in the direction of the derivative. If `is_nonuniform` is false, `d_iv` will take a scalar value of the stepsize between the points used to call the function in `u` and `deriv_iv`. -If `is_nonuniform` is true, the scheme must be able to accept `d_iv` as a vector of stepsizes between the points used to call the function in `u` and `deriv_iv`, therefore of length `length(u)-1`. If it is used on a uniform grid, the stepsizes will be the same for all points. +If `is_nonuniform` is true, the scheme must be able to accept `d_iv` as a vector of stepsizes between the points used to call the function in `u` and `deriv_iv`, therefore of length `length(u)-1`. A method should also be defined for the case where `d_iv` is a scalar, in which case the stepsizes are assumed to be uniform. """ -struct FunctionalScheme{F, V1, V2} <: AbstractScheme +struct FunctionalScheme{F, V1, V2, V3} <: AbstractScheme """ The function that defines the scheme on the interior. """ @@ -75,4 +60,37 @@ struct FunctionalScheme{F, V1, V2} <: AbstractScheme Whether this scheme takes grid steps as input for the nonuniform case. """ is_nonuniform::Bool + """ + parameters for the scheme, should be a vector of numbers. + """ + ps::V3 + """ + The name of the scheme. + """ + name::String +end + +function FunctionalScheme{ips, bps}(interior, lower, upper, is_nonuniform = false, ps = []; name = "FunctionalScheme") where {ips, bps} + + FunctionalScheme{typeof(interior), typeof(lower), + typeof(upper), typeof(ps)}(interior, lower, upper, + ips, bps, is_nonuniform, ps) +end + +function extent(scheme::FunctionalScheme, dorder) + @assert dorder == 1 "Only first order spatial derivatives are implemented for functional schemes." + lower = length(findall(isnothing, scheme.lower)) + upper = length(findall(isnothing, scheme.upper)) + @assert lower == upper "Scheme must have symmetric extent; same number of placeholders in lower and upper boundary functions." + return lower +end + +function lower_extent(scheme::FunctionalScheme, dorder) + @assert dorder == 1 "Only first order spatial derivatives are implemented for functional schemes." + return length(findall(isnothing, scheme.lower)) +end + +function upper_extent(scheme::FunctionalScheme, dorder) + @assert dorder == 1 "Only first order spatial derivatives are implemented for functional schemes." + return length(findall(isnothing, scheme.upper)) end From b63371bf79bbb4ee8b6e79bf41149ff303582dd3 Mon Sep 17 00:00:00 2001 From: xtalax Date: Thu, 2 Mar 2023 20:40:56 +0000 Subject: [PATCH 03/26] fix docs --- docs/pages.jl | 2 +- src/discretization/schemes/WENO/WENO.jl | 2 +- src/interface/scheme_types.jl | 5 ++++- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/docs/pages.jl b/docs/pages.jl index 32c17501e..a7c55f7f1 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -3,7 +3,7 @@ pages = [ "Tutorials" => [ "tutorials/brusselator.md", "tutorials/heat.md", - "tutorials/ivs.md", + "tutorials/params.md", "tutorials/heatss.md", "tutorials/sispde.md", "tutorials/icbc_sampled.md", diff --git a/src/discretization/schemes/WENO/WENO.jl b/src/discretization/schemes/WENO/WENO.jl index ff62f383b..fe6bdaa83 100644 --- a/src/discretization/schemes/WENO/WENO.jl +++ b/src/discretization/schemes/WENO/WENO.jl @@ -58,5 +58,5 @@ end """ function WENOScheme(epsilon = 1e-6) boundry_f = [nothing, nothing] - return FunctionScheme{5, 0}(weno_f, boundary_f, boundary_f, false, [epsilon], name = "WENO") + return FunctionalScheme{5, 0}(weno_f, boundary_f, boundary_f, false, [epsilon], name = "WENO") end diff --git a/src/interface/scheme_types.jl b/src/interface/scheme_types.jl index 2a68e7e23..3ccc670aa 100644 --- a/src/interface/scheme_types.jl +++ b/src/interface/scheme_types.jl @@ -12,7 +12,10 @@ extent(scheme::UpwindScheme, dorder) = dorder + scheme.order - 1 # Functional Schemes """ - # FunctionalScheme +# FunctionalScheme +``` +F = FunctionalScheme{interior_points, boundary_points}(interior, lower, upper, is_nonuniform, parameters; name) +``` A user definable scheme that takes a set of functions as input. The functions define the derivative at the interior, lower boundary, and upper boundary. `lower` and `upper` should be vectors of functions. In general, `upper` and `lower` must be at least `floor(interior_points/2)` long. Where you have no good approximation for a derivative at the boundary, you can use `nothing` as a placeholder. MethodOfLines will then attempt to use an extrapolation here where nessesary. Be warned that this can lead to instability. From e0df0d1aa8027e0a50ee0ac8b72a68172afd848c Mon Sep 17 00:00:00 2001 From: xtalax Date: Thu, 2 Mar 2023 20:42:06 +0000 Subject: [PATCH 04/26] re add docstring --- src/discretization/schemes/WENO/WENO.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/discretization/schemes/WENO/WENO.jl b/src/discretization/schemes/WENO/WENO.jl index fe6bdaa83..08cd526fb 100644 --- a/src/discretization/schemes/WENO/WENO.jl +++ b/src/discretization/schemes/WENO/WENO.jl @@ -1,3 +1,8 @@ +""" +Implements the WENO scheme of Jiang and Shu. +Specified in https://repository.library.brown.edu/studio/item/bdr:297524/PDF/ (Page 8-9) +Implementation *heavily* inspired by https://github.com/ranocha/HyperbolicDiffEq.jl/blob/84c2d882e0c8956457c7d662bf7f18e3c27cfa3d/src/finite_volumes/weno_jiang_shu.jl by H. Ranocha. +""" function weno_f(u, p, t, x, dx) ε = p[1] From 3650241ed0794784deb95298e349eefd958b7a2e Mon Sep 17 00:00:00 2001 From: xtalax Date: Thu, 2 Mar 2023 21:14:09 +0000 Subject: [PATCH 05/26] typo --- src/discretization/schemes/WENO/WENO.jl | 2 +- src/interface/scheme_types.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/discretization/schemes/WENO/WENO.jl b/src/discretization/schemes/WENO/WENO.jl index 08cd526fb..990cda6a4 100644 --- a/src/discretization/schemes/WENO/WENO.jl +++ b/src/discretization/schemes/WENO/WENO.jl @@ -62,6 +62,6 @@ end - `epsilon`: A quantity used to prevent vanishing denominators in the scheme, defaults to `1e-6`. More sensetive problems will benefit from a smaller value. It is defined as a functional scheme. """ function WENOScheme(epsilon = 1e-6) - boundry_f = [nothing, nothing] + boundary_f = [nothing, nothing] return FunctionalScheme{5, 0}(weno_f, boundary_f, boundary_f, false, [epsilon], name = "WENO") end diff --git a/src/interface/scheme_types.jl b/src/interface/scheme_types.jl index 3ccc670aa..52981d3bd 100644 --- a/src/interface/scheme_types.jl +++ b/src/interface/scheme_types.jl @@ -22,7 +22,7 @@ A user definable scheme that takes a set of functions as input. The functions de The boundary functions define the derivative at their index in the function vector, numbering from the boundary. For example, if `boundary_points = 3`, the first function in the vector will define the derivative at the boundary, the second at the boundary plus one step, and the third at the boundary plus two steps. -The functions making up the taking the following inputs: +The functions making up the scheme take the following inputs: Functions must be of the form `f(u, p, t, deriv_iv, d_iv)`. From c8d7312b5fb0f5ec5a764fbc8af242ba41464926 Mon Sep 17 00:00:00 2001 From: xtalax Date: Thu, 2 Mar 2023 21:17:28 +0000 Subject: [PATCH 06/26] properly pass x --- .../schemes/function_scheme/function_scheme.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/discretization/schemes/function_scheme/function_scheme.jl b/src/discretization/schemes/function_scheme/function_scheme.jl index d048ac9b9..aa81cbb9f 100644 --- a/src/discretization/schemes/function_scheme/function_scheme.jl +++ b/src/discretization/schemes/function_scheme/function_scheme.jl @@ -35,18 +35,18 @@ function function_scheme(F::FunctionalScheme{is_nonuniform}, II, s, bs, jx, u, u # Tap points of the stencil, this uses boundary_point_count as this is equal to half the stencil size, which is what we want. u_disc = ufunc(u, Itap, x) ps = params(s) + t = s.time + discx = map(I -> s.grid[x][I[j]], Itap) dx = s.dxs[x] if F.is_nonuniform if dx isa AbstractVector - dx = map(I -> dx[I[j]], Itap)# + dx = map(I -> dx[I[j]], Itap[1:end-1])# end elseif dx isa AbstractVector error("Scheme $(F.name) not implemented for nonuniform dxs.") end - t = s.time - - return f(u_disc, ps, t, x, dx) + return f(u_disc, ps, t, discx, dx) end @inline function generate_advection_rules(F::FunctionalScheme, II::CartesianIndex, s::DiscreteSpace, depvars, derivweights::DifferentialDiscretizer, bcmap, indexmap, terms) From fa9ea0f7f05398bb648a3a1d85f4e29e9395e552 Mon Sep 17 00:00:00 2001 From: xtalax Date: Fri, 3 Mar 2023 15:36:02 +0000 Subject: [PATCH 07/26] assert odd pts --- src/discretization/schemes/function_scheme/function_scheme.jl | 2 +- src/interface/scheme_types.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/discretization/schemes/function_scheme/function_scheme.jl b/src/discretization/schemes/function_scheme/function_scheme.jl index aa81cbb9f..95b04686b 100644 --- a/src/discretization/schemes/function_scheme/function_scheme.jl +++ b/src/discretization/schemes/function_scheme/function_scheme.jl @@ -24,7 +24,7 @@ function get_f_and_taps(F::FunctionalScheme, II, s, bs, jx, u) return f, Itap end -function function_scheme(F::FunctionalScheme{is_nonuniform}, II, s, bs, jx, u, ufunc) where {is_nonuniform<:Val{false}} +function function_scheme(F::FunctionalScheme, II, s, bs, jx, u, ufunc) j, x = jx ndims(u, s) == 0 && return 0 diff --git a/src/interface/scheme_types.jl b/src/interface/scheme_types.jl index 52981d3bd..6d40cdc87 100644 --- a/src/interface/scheme_types.jl +++ b/src/interface/scheme_types.jl @@ -74,7 +74,7 @@ struct FunctionalScheme{F, V1, V2, V3} <: AbstractScheme end function FunctionalScheme{ips, bps}(interior, lower, upper, is_nonuniform = false, ps = []; name = "FunctionalScheme") where {ips, bps} - + @assert ips % 2 == 1 "interior_points must be odd." FunctionalScheme{typeof(interior), typeof(lower), typeof(upper), typeof(ps)}(interior, lower, upper, ips, bps, is_nonuniform, ps) From dbe34b874f601a1fd33235cdf517b626cd084578 Mon Sep 17 00:00:00 2001 From: xtalax Date: Fri, 3 Mar 2023 16:43:52 +0000 Subject: [PATCH 08/26] passing --- docs/src/advection_schemes.md | 6 ++++-- src/MOL_discretization.jl | 4 ++-- src/MOL_utils.jl | 8 ++++++- src/MethodOfLines.jl | 1 + src/discretization/discretize_vars.jl | 6 +++--- .../centered_difference.jl | 2 +- .../function_scheme/function_scheme.jl | 6 +++--- .../integral_expansion/integral_expansion.jl | 4 ++-- .../nonlinear_laplacian.jl | 4 ++-- .../upwind_difference/upwind_difference.jl | 2 +- src/interface/scheme_types.jl | 4 +++- test/pde_systems/MOL_Mixed_Deriv.jl | 6 +++--- test/runtests.jl | 21 +++++++++---------- 13 files changed, 42 insertions(+), 32 deletions(-) diff --git a/docs/src/advection_schemes.md b/docs/src/advection_schemes.md index 97704e673..264cc965e 100644 --- a/docs/src/advection_schemes.md +++ b/docs/src/advection_schemes.md @@ -5,7 +5,7 @@ Used as a keyword argument `advection_scheme` to `MOLFiniteDifference`. ```julia UpwindScheme(approx_order = 1) ``` -Changes the direction of the stencil based on the sign of the coefficient of the odd order derivative to be discretized. Scheme order can be increased by changing the `approx_order` keyword argument. For more information, see [Wikipedia](https://en.wikipedia.org/wiki/Upwind_scheme). +Changes the direction of the stencil based on the sign of the coefficient of the odd order derivative to be discretized. Scheme order can be increased by changing the `approx_order` argument. For more information, see [Wikipedia](https://en.wikipedia.org/wiki/Upwind_scheme). ## WENO Scheme of Jiang and Shu ```julia @@ -21,4 +21,6 @@ Problems with first order derivatives which multiply one another will need to us Supports only first order derivatives, other odd order derivatives are unsupported with this scheme. At present, does not support Nonuniform grids, though this is a planned feature. -Specified on pages 8-9 of [this document](https://repository.library.brown.edu/studio/item/bdr:297524/PDF/) \ No newline at end of file +Specified on pages 8-9 of [this document](https://repository.library.brown.edu/studio/item/bdr:297524/PDF/) + +## FunctionalScheme \ No newline at end of file diff --git a/src/MOL_discretization.jl b/src/MOL_discretization.jl index e884ef99e..4cdfc8418 100644 --- a/src/MOL_discretization.jl +++ b/src/MOL_discretization.jl @@ -4,8 +4,8 @@ function interface_errors(depvars, indvars, discretization) for x in indvars @assert haskey(discretization.dxs, Num(x)) || haskey(discretization.dxs, x) "Variable $x has no step size" end - if !(typeof(discretization.advection_scheme) ∈ [UpwindScheme, WENOScheme]) - throw(ArgumentError("Only `UpwindScheme()` and `WENOScheme()` are supported advection schemes.")) + if !any(s -> discretization.advection_scheme isa s, [UpwindScheme, FunctionalScheme]) + throw(ArgumentError("Only `UpwindScheme()` and `FunctionalScheme()` are supported advection schemes. Got $(typeof(discretization.advection_scheme)).")) end if !(typeof(discretization.disc_strategy) ∈ [ScalarizedDiscretization]) throw(ArgumentError("Only `ScalarizedDiscretization()` are supported discretization strategies.")) diff --git a/src/MOL_utils.jl b/src/MOL_utils.jl index 9cb38d0bd..6aa7851be 100644 --- a/src/MOL_utils.jl +++ b/src/MOL_utils.jl @@ -135,8 +135,14 @@ end function chebyspace(N, dom) interval = dom.domain a, b = DomainSets.infimum(interval), DomainSets.supremum(interval) - x = reverse([(a + b) / 2 + (b - a) / 2 * cos(π * (2k - 1) / (2N)) for k in 1:N]) + x = reverse([(a + b) / 2 + (b - a) / 2 * cospi((2k - 1) / (2N)) for k in 1:N]) x[1] = a x[end] = b return dom.variables => x end + +@inline function sym_dot(a, b) + mapreduce((+), zip(a, b)) do (a_, b_) + a_ * b_ + end +end diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl index 1afac10c8..e282edc8b 100644 --- a/src/MethodOfLines.jl +++ b/src/MethodOfLines.jl @@ -60,6 +60,7 @@ include("system_parsing/pde_system_transformation.jl") include("discretization/interface_boundary.jl") # Schemes +include("discretization/schemes/function_scheme/function_scheme.jl") include("discretization/schemes/centered_difference/centered_difference.jl") include("discretization/schemes/upwind_difference/upwind_difference.jl") include("discretization/schemes/half_offset_centred_difference.jl") diff --git a/src/discretization/discretize_vars.jl b/src/discretization/discretize_vars.jl index 271754df4..5f97f7088 100644 --- a/src/discretization/discretize_vars.jl +++ b/src/discretization/discretize_vars.jl @@ -136,12 +136,12 @@ function DiscreteSpace(vars, discretization::MOLFiniteDifference{G}) where {G} end if t === nothing uaxes = collect(axes(grid[x])[1] for x in arguments(u)) - u => collect(first(@variables $sym[uaxes...])) + u => unwrap.(collect(first(@variables $sym[uaxes...]))) elseif isequal(SymbolicUtils.arguments(u), [t]) - u => fill(u, ()) #Create a 0-dimensional array + u => fill(safe_unwrap(u), ()) #Create a 0-dimensional array else uaxes = collect(axes(grid[x])[1] for x in remove(arguments(u), t)) - u => collect(first(@variables $sym(t)[uaxes...])) + u => unwrap.(collect(first(@variables $sym(t)[uaxes...]))) end end diff --git a/src/discretization/schemes/centered_difference/centered_difference.jl b/src/discretization/schemes/centered_difference/centered_difference.jl index 5869e88b4..f6be5e75a 100644 --- a/src/discretization/schemes/centered_difference/centered_difference.jl +++ b/src/discretization/schemes/centered_difference/centered_difference.jl @@ -23,7 +23,7 @@ function central_difference(D::DerivativeOperator{T,N,Wind,DX}, II, s, bs, jx, u Itap = [bwrap(II + i * I1, bs, s, jx) for i in half_range(D.stencil_length)] end # Tap points of the stencil, this uses boundary_point_count as this is equal to half the stencil size, which is what we want. - return recursive_unwrap(dot(weights, ufunc(u, Itap, x))) + return sym_dot(weights, ufunc(u, Itap, x)) end function central_difference(D::DerivativeOperator{T,N,Wind,DX}, II, s, bs, jx, u, ufunc) where {T,N,Wind,DX<:AbstractVector} diff --git a/src/discretization/schemes/function_scheme/function_scheme.jl b/src/discretization/schemes/function_scheme/function_scheme.jl index 95b04686b..e91430663 100644 --- a/src/discretization/schemes/function_scheme/function_scheme.jl +++ b/src/discretization/schemes/function_scheme/function_scheme.jl @@ -5,8 +5,8 @@ function get_f_and_taps(F::FunctionalScheme, II, s, bs, jx, u) # offset is important due to boundary proximity haslower, hasupper = haslowerupper(bs, x) - lower_point_count = length(f.lower) - upper_point_count = length(f.upper) + lower_point_count = length(F.lower) + upper_point_count = length(F.upper) if (II[j] <= lower_point_count) & !haslower f = F.lower[II[j]] @@ -34,7 +34,7 @@ function function_scheme(F::FunctionalScheme, II, s, bs, jx, u, ufunc) end # Tap points of the stencil, this uses boundary_point_count as this is equal to half the stencil size, which is what we want. u_disc = ufunc(u, Itap, x) - ps = params(s) + ps = vcat(F.ps, params(s)) t = s.time discx = map(I -> s.grid[x][I[j]], Itap) dx = s.dxs[x] diff --git a/src/discretization/schemes/integral_expansion/integral_expansion.jl b/src/discretization/schemes/integral_expansion/integral_expansion.jl index f7bccf73b..b6c60ec68 100644 --- a/src/discretization/schemes/integral_expansion/integral_expansion.jl +++ b/src/discretization/schemes/integral_expansion/integral_expansion.jl @@ -10,7 +10,7 @@ function _euler_integral(II, s, jx, u, ufunc, dx::Number) #where {T,N,Wind,DX<:N Itap = [II - I1, II] weights = [dx / 2, dx / 2] - return dot(weights, ufunc(u, Itap, x)) + _euler_integral(II - I1, s, jx, u, ufunc, dx) + return sym_dot(weights, ufunc(u, Itap, x)) + _euler_integral(II - I1, s, jx, u, ufunc, dx) end # Nonuniform dx @@ -25,7 +25,7 @@ function _euler_integral(II, s, jx, u, ufunc, dx::AbstractVector) #where {T,N,Wi Itap = [II - I1, II] weights = fill(dx[II[j]-1] / 2, 2) - return dot(weights, ufunc(u, Itap, x)) + _euler_integral(II - I1, s, jx, u, ufunc, dx) + return sym_dot(weights, ufunc(u, Itap, x)) + _euler_integral(II - I1, s, jx, u, ufunc, dx) end function euler_integral(II, s, jx, u, ufunc) diff --git a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl index 78d0db9dd..af56d02e6 100644 --- a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl +++ b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl @@ -51,14 +51,14 @@ function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace, interp_weights_and_stencil = [get_half_offset_weights_and_stencil(inner_interpolater, I, s, bs, u, jx) for I in outerstencil] # map variables to symbolically inerpolated/extrapolated expressions - map_vars_to_interpolated(stencil, weights) = [v => dot(weights, s.discvars[v][interface_wrap(stencil)]) for v in depvars] + map_vars_to_interpolated(stencil, weights) = [v => sym_dot(weights, s.discvars[v][interface_wrap(stencil)]) for v in depvars] # Map parameters to interpolated values. Using simplistic extrapolation/interpolation for now as grids are uniform #TODO: make this more efficient map_ivs_to_interpolated(stencil, weights) = safe_vcat([x => dot(weights, getindex.((s.grid[x],), getindex.(interface_wrap(stencil), (j,))))], [s.x̄[k] => s.grid[s.x̄[k]][II[k]] for k in setdiff(1:N, [j])]) # Take the inner finite difference - inner_difference = [dot(inner_weights, s.discvars[u][interface_wrap(inner_stencil)]) for (inner_weights, inner_stencil) in inner_deriv_weights_and_stencil] + inner_difference = [sym_dot(inner_weights, s.discvars[u][interface_wrap(inner_stencil)]) for (inner_weights, inner_stencil) in inner_deriv_weights_and_stencil] # Symbolically interpolate the multiplying expression diff --git a/src/discretization/schemes/upwind_difference/upwind_difference.jl b/src/discretization/schemes/upwind_difference/upwind_difference.jl index 60071e5cb..92626876c 100644 --- a/src/discretization/schemes/upwind_difference/upwind_difference.jl +++ b/src/discretization/schemes/upwind_difference/upwind_difference.jl @@ -65,7 +65,7 @@ function upwind_difference(d::Int, II::CartesianIndex, s::DiscreteSpace, bs, der #@show D.stencil_coefs, D.stencil_length, D.boundary_stencil_length, D.boundary_point_count # unit index in direction of the derivative weights, Itap = _upwind_difference(D, II, s, bs, ispositive, u, jx) - return recursive_unwrap(dot(weights, ufunc(u, Itap, x))) + return sym_dot(weights, ufunc(u, Itap, x)) end function upwind_difference(expr, d::Int, II::CartesianIndex, s::DiscreteSpace, bs, depvars, derivweights, (j, x), u, central_ufunc, indexmap) diff --git a/src/interface/scheme_types.jl b/src/interface/scheme_types.jl index 6d40cdc87..e7efbcaf9 100644 --- a/src/interface/scheme_types.jl +++ b/src/interface/scheme_types.jl @@ -26,6 +26,8 @@ The functions making up the scheme take the following inputs: Functions must be of the form `f(u, p, t, deriv_iv, d_iv)`. +Where the function would branch on a value of u, p, t, or d_iv, use IfElse.ifelse instead of standard conditionals. + For the interior, `u` takes a vector of dependent variable values in the direction of the derivative of length `interior_points`. `interior_points` must be odd, as this function defines the derivative at the center of the input points. @@ -77,7 +79,7 @@ function FunctionalScheme{ips, bps}(interior, lower, upper, is_nonuniform = fals @assert ips % 2 == 1 "interior_points must be odd." FunctionalScheme{typeof(interior), typeof(lower), typeof(upper), typeof(ps)}(interior, lower, upper, - ips, bps, is_nonuniform, ps) + ips, bps, is_nonuniform, ps, name) end function extent(scheme::FunctionalScheme, dorder) diff --git a/test/pde_systems/MOL_Mixed_Deriv.jl b/test/pde_systems/MOL_Mixed_Deriv.jl index f37dc4e09..e923d7db9 100644 --- a/test/pde_systems/MOL_Mixed_Deriv.jl +++ b/test/pde_systems/MOL_Mixed_Deriv.jl @@ -10,7 +10,7 @@ using ModelingToolkit: Differential Dxx = Differential(x)^2 Dtx = Differential(t)*Differential(x) - asf(t, x) = sinpi(2*(t + (1+sqrt(5))*x/2)) + cospi(2*(t + (1-sqrt(5))*x/2)) + assf(t, x) = sinpi(2*(t + (1+sqrt(5))*x/2)) + cospi(2*(t + (1-sqrt(5))*x/2)) aDtf(t, x) = 2pi*cospi(2*(t + (1+sqrt(5))*x/2)) - 2pi*sinpi(2*(t + (1-sqrt(5))*x/2)) # Where asf(t, x) ~ 0, NonlinearSolved xmin = -0.1118033987645 @@ -18,7 +18,7 @@ using ModelingToolkit: Differential eq = [Dtt(u(t, x)) + Dtx(u(t, x)) - Dxx(u(t, x)) ~ Dxx(x)] - bcs = [u(1e-9, x) ~ asf(1e-9, x), + bcs = [u(1e-9, x) ~ assf(1e-9, x), Dt(u(1e-9, x)) ~ aDtf(1e-9, x), u(t, xmin) ~ 0, u(t, xmax) ~ 0] @@ -39,7 +39,7 @@ using ModelingToolkit: Differential tdisc = sol[t] usol = sol[u(t,x)] - asol = [asf(t, x) for t in tdisc, x in xdisc] + asol = [assf(t, x) for t in tdisc, x in xdisc] @test_broken usol ≈ asol atol = 1e-3 end diff --git a/test/runtests.jl b/test/runtests.jl index 65c0cab7e..a730b1ad7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,11 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") # Start Test Script @time begin + if GROUP == "All" || GROUP == "Convection_WENO" + @time @safetestset "MOLFiniteDifference Interface: Linear Convection, WENO Scheme." begin + include("pde_systems/MOL_1D_Linear_Convection_WENO.jl") + end + end if GROUP == "All" || GROUP == "Diffusion_NU" @time @safetestset "MOLFiniteDifference Interface: 1D Linear Diffusion, Non-Uniform" begin @@ -100,12 +105,6 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") end end - if GROUP == "All" || GROUP == "Convection_WENO" - @time @safetestset "MOLFiniteDifference Interface: Linear Convection, WENO Scheme." begin - include("pde_systems/MOL_1D_Linear_Convection_WENO.jl") - end - end - if GROUP == "All" || GROUP == "Burgers" @time @safetestset "MOLFiniteDifference Interface: 2D Burger's Equation" begin include("pde_systems/burgers_eq.jl") @@ -124,9 +123,9 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") end end - # if GROUP == "All" || GROUP == "Mixed_Derivatives" - # @time @safetestset "MOLFiniteDifference Interface: Mixed Derivatives" begin - # include("pde_systems/MOL_Mixed_Deriv.jl") - # end - # end + if GROUP == "All" || GROUP == "Mixed_Derivatives" + @time @safetestset "MOLFiniteDifference Interface: Mixed Derivatives" begin + include("pde_systems/MOL_Mixed_Deriv.jl") + end + end end From 6d32becaa7f6836adab9757a3729e5619c5b726c Mon Sep 17 00:00:00 2001 From: xtalax Date: Mon, 6 Mar 2023 15:34:32 +0000 Subject: [PATCH 09/26] more sym_dot --- .../schemes/centered_difference/centered_difference.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/discretization/schemes/centered_difference/centered_difference.jl b/src/discretization/schemes/centered_difference/centered_difference.jl index f6be5e75a..b47d6dc62 100644 --- a/src/discretization/schemes/centered_difference/centered_difference.jl +++ b/src/discretization/schemes/centered_difference/centered_difference.jl @@ -48,7 +48,7 @@ function central_difference(D::DerivativeOperator{T,N,Wind,DX}, II, s, bs, jx, u end # Tap points of the stencil, this uses boundary_point_count as this is equal to half the stencil size, which is what we want. - return recursive_unwrap(dot(weights, ufunc(u, Itap, x))) + return sym_dot(weights, ufunc(u, Itap, x)) end """ From f60258be288abab21cc90126cd60c0183a99e9e7 Mon Sep 17 00:00:00 2001 From: xtalax Date: Mon, 6 Mar 2023 15:38:29 +0000 Subject: [PATCH 10/26] new retcode --- test/pde_systems/MOL_1D_HigherOrder.jl | 2 +- test/pde_systems/MOL_1D_NonLinear_Diffusion.jl | 14 +++++++------- .../MOL_1D_NonLinear_Diffusion_NonUniform.jl | 14 +++++++------- test/pde_systems/MOL_Mixed_Deriv.jl | 2 +- test/pde_systems/MOLtest2.jl | 2 +- 5 files changed, 17 insertions(+), 17 deletions(-) diff --git a/test/pde_systems/MOL_1D_HigherOrder.jl b/test/pde_systems/MOL_1D_HigherOrder.jl index 83fb6f479..a21aeab5d 100644 --- a/test/pde_systems/MOL_1D_HigherOrder.jl +++ b/test/pde_systems/MOL_1D_HigherOrder.jl @@ -114,7 +114,7 @@ end sol = solve(prob,FBDF(),saveat=0.1) - @test sol.retcode == :Success + @test sol.retcode == SciMLBase.ReturnCode.Success xs = sol[x] ts = sol[t] diff --git a/test/pde_systems/MOL_1D_NonLinear_Diffusion.jl b/test/pde_systems/MOL_1D_NonLinear_Diffusion.jl index 2c569372a..bf925fae7 100644 --- a/test/pde_systems/MOL_1D_NonLinear_Diffusion.jl +++ b/test/pde_systems/MOL_1D_NonLinear_Diffusion.jl @@ -46,7 +46,7 @@ using ModelingToolkit: Differential # Solution of the ODE system using OrdinaryDiffEq sol = solve(prob, Rosenbrock32()) - @test sol.retcode == :Success + @test sol.retcode == SciMLBase.ReturnCode.Success # Test against exact solution x_disc = sol[x] @@ -161,7 +161,7 @@ end # Solution of the ODE system using OrdinaryDiffEq sol = solve(prob, Rosenbrock32()) - @test sol.retcode == :Success + @test sol.retcode == SciMLBase.ReturnCode.Success # Test against exact solution x_disc = sol[x] @@ -222,7 +222,7 @@ end using OrdinaryDiffEq sol = solve(prob, Rosenbrock32()) - @test sol.retcode == :Success + @test sol.retcode == SciMLBase.ReturnCode.Success # Test against exact solution x_disc = sol[x] @@ -337,7 +337,7 @@ end # Solution of the ODE system using OrdinaryDiffEq sol = solve(prob, Rosenbrock32()) # TODO: check warnings - @test sol.retcode == :Success + @test sol.retcode == SciMLBase.ReturnCode.Success # Test against exact solution x_disc = sol[x] @@ -394,7 +394,7 @@ end # Solution of the ODE system using OrdinaryDiffEq sol = solve(prob, Rosenbrock32()) # TODO: check warnings - @test sol.retcode == :Success + @test sol.retcode == SciMLBase.ReturnCode.Success # Test against exact solution @@ -453,7 +453,7 @@ end # Solution of the ODE system using OrdinaryDiffEq sol = solve(prob, Rosenbrock32()) - @test_broken sol.retcode == :Success + @test_broken sol.retcode == SciMLBase.ReturnCode.Success fail # don't run the rest of the test # Test against exact solution x_disc = sol[x] @@ -510,7 +510,7 @@ end # Solution of the ODE system using OrdinaryDiffEq sol = solve(prob, Rosenbrock32()) - @test_broken sol.retcode == :Success + @test_broken sol.retcode == SciMLBase.ReturnCode.Success fail #don't run later tests # Test against exact solution x_disc = sol[x] diff --git a/test/pde_systems/MOL_1D_NonLinear_Diffusion_NonUniform.jl b/test/pde_systems/MOL_1D_NonLinear_Diffusion_NonUniform.jl index db4959296..1f75c76ce 100644 --- a/test/pde_systems/MOL_1D_NonLinear_Diffusion_NonUniform.jl +++ b/test/pde_systems/MOL_1D_NonLinear_Diffusion_NonUniform.jl @@ -52,7 +52,7 @@ using StableRNGs using OrdinaryDiffEq sol = solve(prob, Rosenbrock32()) - @test sol.retcode == :Success + @test sol.retcode == SciMLBase.ReturnCode.Success # Test against exact solution x_disc = sol[x] @@ -173,7 +173,7 @@ end using OrdinaryDiffEq sol = solve(prob, Rosenbrock32()) - @test sol.retcode == :Success + @test sol.retcode == SciMLBase.ReturnCode.Success # Test against exact solution x_disc = sol[x] @@ -236,7 +236,7 @@ end using OrdinaryDiffEq sol = solve(prob, Rosenbrock32()) - @test sol.retcode == :Success + @test sol.retcode == SciMLBase.ReturnCode.Success # Test against exact solution x_disc = sol[x] @@ -353,7 +353,7 @@ end # Solution of the ODE system using OrdinaryDiffEq sol = solve(prob, Rosenbrock32()) - @test sol.retcode == :Success + @test sol.retcode == SciMLBase.ReturnCode.Success # Test against exact solution x_disc = sol[x] @@ -414,7 +414,7 @@ end using OrdinaryDiffEq sol = solve(prob, Rosenbrock32()) - @test sol.retcode == :Success + @test sol.retcode == SciMLBase.ReturnCode.Success # Test against exact solution x_disc = sol[x] @@ -476,7 +476,7 @@ end using OrdinaryDiffEq sol = solve(prob, Rosenbrock32()) - @test_broken sol.retcode == :Success + @test_broken sol.retcode == SciMLBase.ReturnCode.Success fail # Test against exact solution x_disc = sol[x] @@ -538,7 +538,7 @@ end using OrdinaryDiffEq sol = solve(prob, Rosenbrock32()) - @test_broken sol.retcode == :Success + @test_broken sol.retcode == SciMLBase.ReturnCode.Success fail # Test against exact solution x_disc = sol[x] diff --git a/test/pde_systems/MOL_Mixed_Deriv.jl b/test/pde_systems/MOL_Mixed_Deriv.jl index e923d7db9..50b90965a 100644 --- a/test/pde_systems/MOL_Mixed_Deriv.jl +++ b/test/pde_systems/MOL_Mixed_Deriv.jl @@ -67,5 +67,5 @@ end prob = discretize(pdesys, discretization, advection_scheme = WENOScheme()) sol = solve(prob, FBDF(), saveat = 0.1); - @test sol.retcode == :Success + @test sol.retcode == SciMLBase.ReturnCode.Successend end diff --git a/test/pde_systems/MOLtest2.jl b/test/pde_systems/MOLtest2.jl index e87496fdb..cef11ec27 100644 --- a/test/pde_systems/MOLtest2.jl +++ b/test/pde_systems/MOLtest2.jl @@ -156,7 +156,7 @@ end solc = hcat(solc1[:, :], solc2[:, 2:end]) - @test sol.retcode == :Success + @test sol.retcode == SciMLBase.ReturnCode.Success end @testset "Another boundaries appearing in equations case" begin From bc0da37d7b74ccbdc134305bf3bed2400e666e70 Mon Sep 17 00:00:00 2001 From: xtalax Date: Mon, 6 Mar 2023 17:43:07 +0000 Subject: [PATCH 11/26] more sym_dot --- .../schemes/nonlinear_laplacian/nonlinear_laplacian.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl index af56d02e6..7360c7acd 100644 --- a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl +++ b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl @@ -68,7 +68,7 @@ function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace, end # multiply the inner finite difference by the interpolated expression, and finally take the outer finite difference - return dot(outerweights, inner_difference .* interpolated_expr) + return sym_dot(outerweights, inner_difference .* interpolated_expr) end @inline function generate_nonlinlap_rules(II::CartesianIndex, s::DiscreteSpace, depvars, derivweights::DifferentialDiscretizer, bcmap, indexmap, terms) From 122e2f9e0ae9ee0b68c6947f7e7e0976ca720541 Mon Sep 17 00:00:00 2001 From: xtalax Date: Mon, 6 Mar 2023 18:27:49 +0000 Subject: [PATCH 12/26] ho dot --- src/discretization/schemes/half_offset_centred_difference.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/discretization/schemes/half_offset_centred_difference.jl b/src/discretization/schemes/half_offset_centred_difference.jl index 2ae1a4357..8b6a00a8f 100644 --- a/src/discretization/schemes/half_offset_centred_difference.jl +++ b/src/discretization/schemes/half_offset_centred_difference.jl @@ -61,5 +61,5 @@ function half_offset_centered_difference(D, II, s, bs, jx, u, ufunc) ndims(u, s) == 0 && return Num(0) j, x = jx weights, Itap = get_half_offset_weights_and_stencil(D, II, s, bs, u, jx) - return recursive_unwrap(dot(weights, ufunc(u, Itap, x))) + return sym_dot(weights, ufunc(u, Itap, x)) end From 38210d5a3dc0bddc4c2c80f56aca5ab1064accb1 Mon Sep 17 00:00:00 2001 From: xtalax Date: Fri, 10 Mar 2023 11:56:45 +0000 Subject: [PATCH 13/26] fix extrap x interface --- src/discretization/generate_bc_eqs.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/discretization/generate_bc_eqs.jl b/src/discretization/generate_bc_eqs.jl index f94d508ef..f72b00806 100644 --- a/src/discretization/generate_bc_eqs.jl +++ b/src/discretization/generate_bc_eqs.jl @@ -167,7 +167,12 @@ function generate_extrap_eqs!(eqs, pde, u, s, derivweights, interiormap, bcmap) for (j, x) in enumerate(args) ninterp = lowerextents[j] - vlower[j] I1 = unitindex(length(args), j) + bs = bcmap[operation(u)][x] + haslower, hasupper = haslowerupper(bs, x) while ninterp >= vlower[j] + if haslower + break + end for Il in (edge(interiormap, s, u, j, true) .+ (ninterp * I1,)) expr = central_difference(derivweights.boundary[x], Il, s, filter_interfaces(bcmap[operation(u)][x]), (j, x), u, ufunc) push!(eqmap[Il], expr) @@ -176,6 +181,9 @@ function generate_extrap_eqs!(eqs, pde, u, s, derivweights, interiormap, bcmap) end ninterp = upperextents[j] - vupper[j] while ninterp >= vupper[j] + if hasupper + break + end for Iu in (edge(interiormap, s, u, j, false) .- (ninterp * I1,)) expr = central_difference(derivweights.boundary[x], Iu, s, filter_interfaces(bcmap[operation(u)][x]), (j, x), u, ufunc) push!(eqmap[Iu], expr) From cfbc2b0a137f84bf32ece7b80d6b39d9b534b6fb Mon Sep 17 00:00:00 2001 From: xtalax Date: Fri, 10 Mar 2023 11:56:58 +0000 Subject: [PATCH 14/26] new analytic sig --- src/MOL_utils.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/MOL_utils.jl b/src/MOL_utils.jl index 6aa7851be..202e72400 100644 --- a/src/MOL_utils.jl +++ b/src/MOL_utils.jl @@ -84,14 +84,14 @@ function generate_function_from_gridlocs(analyticmap, gridlocs, s) is_t_first = is_t_first_map[uop] _f = analyticmap[opsmap[uop]] if is_t_first - return t -> _f(t, x̄...) + return (p, t) -> _f(p, t, x̄...) else - return t -> _f(x̄..., t) + return (p, t) -> _f(p, x̄..., t) end end f = (u0, p, t) -> map(fs_) do f_ - f_(t) + f_(p, t) end return f From 7229b8a647fe763cb9051094eeeb26a48a524eeb Mon Sep 17 00:00:00 2001 From: xtalax Date: Fri, 10 Mar 2023 11:57:12 +0000 Subject: [PATCH 15/26] update test --- test/components/ODEFunction_test.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/components/ODEFunction_test.jl b/test/components/ODEFunction_test.jl index ecb701126..15516a09b 100644 --- a/test/components/ODEFunction_test.jl +++ b/test/components/ODEFunction_test.jl @@ -1,6 +1,6 @@ using ModelingToolkit, MethodOfLines, DomainSets, OrdinaryDiffEq, Test -@parameters x t +@parameters x t a @variables u(..) Dx = Differential(x) Dt = Differential(t) @@ -9,7 +9,7 @@ x_max = 1.0 t_min = 0.0 t_max = 6.0 -analytic_u(t, x) = x / (t + 1) +analytic_u(p, t, x) = x / (t + p[1]) eq = Dt(u(t, x)) ~ -u(t, x) * Dx(u(t, x)) @@ -20,11 +20,11 @@ bcs = [u(0, x) ~ x, domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)] -@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)]) +@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)], [a => 1.0]) disc = MOLFiniteDifference([x => 30], t, advection_scheme=WENOScheme()) -prob = discretize(pdesys, disc; analytic = [u(t, x) => analytic_u]) +prob = discretize(pdesys, disc; analytic = pdesys.analytic_func) sol = solve(prob, FBDF()) From 4a0b61e018fa9ce65892969d6955338c3c001910 Mon Sep 17 00:00:00 2001 From: xtalax Date: Fri, 10 Mar 2023 12:50:43 +0000 Subject: [PATCH 16/26] fix test --- test/components/ODEFunction_test.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/components/ODEFunction_test.jl b/test/components/ODEFunction_test.jl index 15516a09b..9f5f51030 100644 --- a/test/components/ODEFunction_test.jl +++ b/test/components/ODEFunction_test.jl @@ -14,13 +14,13 @@ analytic_u(p, t, x) = x / (t + p[1]) eq = Dt(u(t, x)) ~ -u(t, x) * Dx(u(t, x)) bcs = [u(0, x) ~ x, - u(t, x_min) ~ analytic_u(t, x_min), - u(t, x_max) ~ analytic_u(t, x_max)] + u(t, x_min) ~ analytic_u(_, t, x_min), + u(t, x_max) ~ analytic_u(_, t, x_max)] domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)] -@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)], [a => 1.0]) +@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)], [a => 1.0], analytic_func = [u(t, x) => analytic_u]) disc = MOLFiniteDifference([x => 30], t, advection_scheme=WENOScheme()) From 246cbf71694f4fec73a56f19275488c838d1ec35 Mon Sep 17 00:00:00 2001 From: xtalax Date: Fri, 10 Mar 2023 12:54:45 +0000 Subject: [PATCH 17/26] fix test better --- test/components/ODEFunction_test.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/components/ODEFunction_test.jl b/test/components/ODEFunction_test.jl index 9f5f51030..b68133f8b 100644 --- a/test/components/ODEFunction_test.jl +++ b/test/components/ODEFunction_test.jl @@ -14,8 +14,8 @@ analytic_u(p, t, x) = x / (t + p[1]) eq = Dt(u(t, x)) ~ -u(t, x) * Dx(u(t, x)) bcs = [u(0, x) ~ x, - u(t, x_min) ~ analytic_u(_, t, x_min), - u(t, x_max) ~ analytic_u(_, t, x_max)] + u(t, x_min) ~ analytic_u([1], t, x_min), + u(t, x_max) ~ analytic_u([1], t, x_max)] domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)] From 1784cb18592950f1b44af30a06d3451a040863ab Mon Sep 17 00:00:00 2001 From: xtalax Date: Fri, 10 Mar 2023 13:03:30 +0000 Subject: [PATCH 18/26] a comment --- src/MOL_utils.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/MOL_utils.jl b/src/MOL_utils.jl index 202e72400..6b482fe30 100644 --- a/src/MOL_utils.jl +++ b/src/MOL_utils.jl @@ -64,6 +64,7 @@ end function get_gridloc(u, s) if isequal(operation(u), getindex) + # Remember arguments of getindex have u(t) first return _get_gridloc(s, arguments(u)...) else return (operation(u), []) From 6331a3b59b3ed43abd35da9a30cd978d3ca551b4 Mon Sep 17 00:00:00 2001 From: xtalax Date: Mon, 13 Mar 2023 14:55:50 +0000 Subject: [PATCH 19/26] views --- .../schemes/function_scheme/function_scheme.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/discretization/schemes/function_scheme/function_scheme.jl b/src/discretization/schemes/function_scheme/function_scheme.jl index e91430663..757289e72 100644 --- a/src/discretization/schemes/function_scheme/function_scheme.jl +++ b/src/discretization/schemes/function_scheme/function_scheme.jl @@ -36,11 +36,12 @@ function function_scheme(F::FunctionalScheme, II, s, bs, jx, u, ufunc) u_disc = ufunc(u, Itap, x) ps = vcat(F.ps, params(s)) t = s.time - discx = map(I -> s.grid[x][I[j]], Itap) + itap = map(I -> I[j], Itap) + discx = @view s.grid[x][itap] dx = s.dxs[x] if F.is_nonuniform if dx isa AbstractVector - dx = map(I -> dx[I[j]], Itap[1:end-1])# + dx = @views dx[itap[1:end-1]] end elseif dx isa AbstractVector error("Scheme $(F.name) not implemented for nonuniform dxs.") From 432baaa683858945b2cb6d8ed9d5ac1ea0e33040 Mon Sep 17 00:00:00 2001 From: xtalax Date: Mon, 13 Mar 2023 15:13:33 +0000 Subject: [PATCH 20/26] update stats --- Project.toml | 2 +- src/interface/solution/timedep.jl | 4 ++-- src/interface/solution/timeindep.jl | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index a50d1058a..acf102f7d 100644 --- a/Project.toml +++ b/Project.toml @@ -29,7 +29,7 @@ Latexify = "0.15" ModelingToolkit = "8" OrdinaryDiffEq = "6" SafeTestsets = "0.0.1" -SciMLBase = "1" +SciMLBase = "1.90" StaticArrays = "1" SymbolicUtils = "1" Symbolics = "5" diff --git a/src/interface/solution/timedep.jl b/src/interface/solution/timedep.jl index 8bad97a41..f0c731ed6 100644 --- a/src/interface/solution/timedep.jl +++ b/src/interface/solution/timedep.jl @@ -50,10 +50,10 @@ function SciMLBase.PDETimeSeriesSolution(sol::SciMLBase.AbstractODESolution{T}, return SciMLBase.PDETimeSeriesSolution{T,length(discretespace.ū),typeof(umap),typeof(metadata), typeof(sol),typeof(sol.errors),typeof(sol.t),typeof(ivgrid), typeof(ivs),typeof(pdesys.dvs),typeof(sol.prob),typeof(sol.alg), - typeof(interp)}(umap, sol, sol.errors, sol.t, ivgrid, ivs, + typeof(interp), typeof(sol.stats)}(umap, sol, sol.errors, sol.t, ivgrid, ivs, pdesys.dvs, metadata, sol.prob, sol.alg, interp, sol.dense, sol.tslocation, - sol.retcode) + sol.retcode, sol.stats) catch e rethrow(e) return sol, e diff --git a/src/interface/solution/timeindep.jl b/src/interface/solution/timeindep.jl index fbd5056c0..3fce94195 100644 --- a/src/interface/solution/timeindep.jl +++ b/src/interface/solution/timeindep.jl @@ -30,7 +30,7 @@ function SciMLBase.PDENoTimeSolution(sol::SciMLBase.NonlinearSolution{T}, metada return SciMLBase.PDENoTimeSolution{T,length(discretespace.ū),typeof(umap),typeof(metadata), typeof(sol),typeof(ivgrid),typeof(ivs),typeof(pdesys.dvs),typeof(sol.prob),typeof(sol.alg), - typeof(interp)}(umap, sol, ivgrid, ivs, + typeof(interp), typeof(sol.stats)}(umap, sol, ivgrid, ivs, pdesys.dvs, metadata, sol.prob, sol.alg, - interp, sol.retcode) + interp, sol.retcode, sol.stats) end From 0b6350dcf7ec12d9121c151b516f8d3da995d305 Mon Sep 17 00:00:00 2001 From: xtalax Date: Mon, 20 Mar 2023 14:40:29 +0000 Subject: [PATCH 21/26] fix odefunc --- src/MOL_discretization.jl | 3 ++- src/MOL_utils.jl | 1 - src/MethodOfLines.jl | 5 +++-- src/system_parsing/variable_map.jl | 7 +++++-- test/components/ODEFunction_test.jl | 2 +- test/pde_systems/MOL_Mixed_Deriv.jl | 2 +- test/runtests.jl | 25 +++++++++++++------------ 7 files changed, 25 insertions(+), 20 deletions(-) diff --git a/src/MOL_discretization.jl b/src/MOL_discretization.jl index 4cdfc8418..1312db733 100644 --- a/src/MOL_discretization.jl +++ b/src/MOL_discretization.jl @@ -147,7 +147,8 @@ function SciMLBase.discretize(pdesys::PDESystem,discretization::MethodOfLines.MO return prob else f = ODEFunction(pdesys, discretization, analytic=analytic, discretization.kwargs..., kwargs...) - return ODEProblem(f, prob.u0, prob.tspan; discretization.kwargs..., kwargs...) + + return ODEProblem(f, prob.u0, prob.tspan, prob.p; discretization.kwargs..., kwargs...) end end end diff --git a/src/MOL_utils.jl b/src/MOL_utils.jl index 6b482fe30..2fcc11deb 100644 --- a/src/MOL_utils.jl +++ b/src/MOL_utils.jl @@ -71,7 +71,6 @@ function get_gridloc(u, s) end end - function generate_function_from_gridlocs(analyticmap, gridlocs, s) is_t_first_map = Dict(map(s.ū) do u operation(u) => (findfirst(x -> isequal(s.time, x), arguments(u)) == 1) diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl index e282edc8b..5b8296401 100644 --- a/src/MethodOfLines.jl +++ b/src/MethodOfLines.jl @@ -3,7 +3,8 @@ using LinearAlgebra using SciMLBase using DiffEqBase using ModelingToolkit -using ModelingToolkit: operation, istree, arguments, variable, get_metadata, get_states +using ModelingToolkit: operation, istree, arguments, variable, get_metadata, get_states, +parameters, defaults, varmap_to_vars using SymbolicUtils, Symbolics using Symbolics: unwrap, solve_for, expand_derivatives, diff2term, setname, rename, similarterm using SymbolicUtils: operation, arguments @@ -80,6 +81,6 @@ include("scalar_discretization.jl") include("MOL_discretization.jl") export MOLFiniteDifference, discretize, symbolic_discretize, ODEFunctionExpr, generate_code, grid_align, edge_align, center_align, get_discrete, chebyspace -export UpwindScheme, WENOScheme +export UpwindScheme, WENOScheme, FunctionalScheme end diff --git a/src/system_parsing/variable_map.jl b/src/system_parsing/variable_map.jl index 4d8bdb573..8e3fd3ed3 100644 --- a/src/system_parsing/variable_map.jl +++ b/src/system_parsing/variable_map.jl @@ -11,6 +11,9 @@ struct VariableMap end function VariableMap(eqs, depvars, domain, time, ps = []) + if ps isa SciMLBase.NullParameters + ps = [] + end time = safe_unwrap(time) ps = map(ps) do p safe_unwrap(p.first) @@ -34,8 +37,8 @@ function VariableMap(eqs, depvars, domain, time, ps = []) return VariableMap(ū, x̄, ps, time, Dict(intervals), Dict(args), depvar_ops, Dict(x̄2dim), Dict(dim2x̄)) end -VariableMap(pdesys::PDESystem, disc::MOLFiniteDifference) = VariableMap(pdesys.eqs, pdesys.dvs, pdesys.domain, disc.time) -VariableMap(pdesys::PDESystem, t) = VariableMap(pdesys.eqs, pdesys.dvs, pdesys.domain, t) +VariableMap(pdesys::PDESystem, disc::MOLFiniteDifference) = VariableMap(pdesys.eqs, pdesys.dvs, pdesys.domain, disc.time, pdesys.ps) +VariableMap(pdesys::PDESystem, t) = VariableMap(pdesys.eqs, pdesys.dvs, pdesys.domain, t, pdesys.ps) function update_varmap!(v, newdv) push!(v.ū, newdv) diff --git a/test/components/ODEFunction_test.jl b/test/components/ODEFunction_test.jl index b68133f8b..b0cf165d8 100644 --- a/test/components/ODEFunction_test.jl +++ b/test/components/ODEFunction_test.jl @@ -11,7 +11,7 @@ t_max = 6.0 analytic_u(p, t, x) = x / (t + p[1]) -eq = Dt(u(t, x)) ~ -u(t, x) * Dx(u(t, x)) +eq = Dt(u(t, x)) ~ -a*u(t, x) * Dx(u(t, x)) bcs = [u(0, x) ~ x, u(t, x_min) ~ analytic_u([1], t, x_min), diff --git a/test/pde_systems/MOL_Mixed_Deriv.jl b/test/pde_systems/MOL_Mixed_Deriv.jl index 50b90965a..dde8d1a29 100644 --- a/test/pde_systems/MOL_Mixed_Deriv.jl +++ b/test/pde_systems/MOL_Mixed_Deriv.jl @@ -67,5 +67,5 @@ end prob = discretize(pdesys, discretization, advection_scheme = WENOScheme()) sol = solve(prob, FBDF(), saveat = 0.1); - @test sol.retcode == SciMLBase.ReturnCode.Successend + @test sol.retcode == SciMLBase.ReturnCode.Success end diff --git a/test/runtests.jl b/test/runtests.jl index a730b1ad7..a78770afd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,18 +7,6 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") # Start Test Script @time begin - if GROUP == "All" || GROUP == "Convection_WENO" - @time @safetestset "MOLFiniteDifference Interface: Linear Convection, WENO Scheme." begin - include("pde_systems/MOL_1D_Linear_Convection_WENO.jl") - end - end - - if GROUP == "All" || GROUP == "Diffusion_NU" - @time @safetestset "MOLFiniteDifference Interface: 1D Linear Diffusion, Non-Uniform" begin - include("pde_systems/MOL_1D_Linear_Diffusion_NonUniform.jl") - end - end - if GROUP == "All" || GROUP == "Components" @time @safetestset "MOLFiniteDifference Utils" begin include("utils_test.jl") @@ -38,6 +26,19 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") #@time @safetestset "Finite Difference Schemes" begin include("components/finite_diff_schemes.jl") end end + if GROUP == "All" || GROUP == "Convection_WENO" + @time @safetestset "MOLFiniteDifference Interface: Linear Convection, WENO Scheme." begin + include("pde_systems/MOL_1D_Linear_Convection_WENO.jl") + end + end + + if GROUP == "All" || GROUP == "Diffusion_NU" + @time @safetestset "MOLFiniteDifference Interface: 1D Linear Diffusion, Non-Uniform" begin + include("pde_systems/MOL_1D_Linear_Diffusion_NonUniform.jl") + end + end + + if GROUP == "All" || GROUP == "Integrals" @time @safetestset "MOLFiniteDifference Interface: Integrals" begin include("pde_systems/MOL_1D_Integration.jl") From 4d0391127d79943d35b7afa68d56debc96f35505 Mon Sep 17 00:00:00 2001 From: xtalax Date: Tue, 21 Mar 2023 15:08:08 +0000 Subject: [PATCH 22/26] fix kwarg --- src/discretization/generate_bc_eqs.jl | 2 +- src/scalar_discretization.jl | 10 ++++++++-- src/system_parsing/interior_map.jl | 2 +- 3 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/discretization/generate_bc_eqs.jl b/src/discretization/generate_bc_eqs.jl index f72b00806..43ff62dd8 100644 --- a/src/discretization/generate_bc_eqs.jl +++ b/src/discretization/generate_bc_eqs.jl @@ -209,7 +209,7 @@ end @inline function generate_corner_eqs!(bceqs, s, interiormap, N, u) interior = interiormap.I[interiormap.pde[u]] - ndims(u, s) == 0 && return + N == 0 || ndims(u, s) == 0 && return sd(i, j) = selectdim(interior, j, i) domain = setdiff(s.Igrid[u], interior) II1 = unitindices(N) diff --git a/src/scalar_discretization.jl b/src/scalar_discretization.jl index 2ecb1f840..04dbbf56b 100644 --- a/src/scalar_discretization.jl +++ b/src/scalar_discretization.jl @@ -47,18 +47,24 @@ function generate_system(alleqs, bceqs, ics, discvars, defaults, ps, tspan, meta alleqs = vcat(alleqs, unique(bceqs)) alldepvarsdisc = vec(reduce(vcat, vec(unique(reduce(vcat, vec.(values(discvars))))))) # Finalize + if haskey(metadata.disc.kwargs, :checks) + checks = metadata.disc.kwargs[:checks] + else + checks = true + end + end try if t === nothing # At the time of writing, NonlinearProblems require that the system of equations be in this form: # 0 ~ ... # Thus, before creating a NonlinearSystem we normalize the equations s.t. the lhs is zero. eqs = map(eq -> 0 ~ eq.rhs - eq.lhs, alleqs) - sys = NonlinearSystem(eqs, alldepvarsdisc, ps, defaults=defaults, name=name, metadata=metadata) + sys = NonlinearSystem(eqs, alldepvarsdisc, ps, defaults=defaults, name=name, metadata=metadata, checks = checks) return sys, nothing else # * In the end we have reduced the problem to a system of equations in terms of Dt that can be solved by an ODE solver. - sys = ODESystem(alleqs, t, alldepvarsdisc, ps, defaults=defaults, name=name, metadata=metadata) + sys = ODESystem(alleqs, t, alldepvarsdisc, ps, defaults=defaults, name=name, metadata=metadata, checks = checks) return sys, tspan end catch e diff --git a/src/system_parsing/interior_map.jl b/src/system_parsing/interior_map.jl index 1fdfccf7a..49790f600 100644 --- a/src/system_parsing/interior_map.jl +++ b/src/system_parsing/interior_map.jl @@ -59,7 +59,7 @@ end function generate_interior(lower, upper, u, s, disc::MOLFiniteDifference{G,D}) where {G, D<:ScalarizedDiscretization} args = remove(arguments(u), s.time) - ret = s.Igrid[u][[(1+lower[x2i(s, u, x)]:length(s.grid[x])-upper[x2i(s, u, x)]) for x in args]...] + ret = collect(s.Igrid[u][[(1+lower[x2i(s, u, x)]:length(s.grid[x])-upper[x2i(s, u, x)]) for x in args]...]) return ret end From e5078267748f3af3be92409fe3f7e05b92dc6825 Mon Sep 17 00:00:00 2001 From: xtalax Date: Tue, 21 Mar 2023 15:30:03 +0000 Subject: [PATCH 23/26] fix docs --- Project.toml | 2 +- docs/src/tutorials/PIDE.md | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index acf102f7d..3b8b5abd4 100644 --- a/Project.toml +++ b/Project.toml @@ -29,7 +29,7 @@ Latexify = "0.15" ModelingToolkit = "8" OrdinaryDiffEq = "6" SafeTestsets = "0.0.1" -SciMLBase = "1.90" +SciMLBase = "1.90, 1.91" StaticArrays = "1" SymbolicUtils = "1" Symbolics = "5" diff --git a/docs/src/tutorials/PIDE.md b/docs/src/tutorials/PIDE.md index 18c946da4..1dcd16b3b 100644 --- a/docs/src/tutorials/PIDE.md +++ b/docs/src/tutorials/PIDE.md @@ -56,7 +56,7 @@ To have an integral over the whole domain, be sure to wrap the integral in an au Due to a limitation, the whole domain integral needs to have the same arguments as the integrand, but is constant in x. To use it in an equation one dimension lower, use a boundary value like integral(t, 0.0) ```@example integrals2 -using MethodOfLines, ModelingToolkit, DomainSets +using MethodOfLines, ModelingToolkit, DomainSets, OrdinaryDiffEq, Plots @parameters t, x @variables integrand(..) integral(..) @@ -90,5 +90,6 @@ xdisc = sol[x] tdisc = sol[t] integralsol = sol[integral(t, x)] -exact = [asf(t_) for t_ in tdisc, x_ in xdisc] + +heatmap(integralsol) ``` \ No newline at end of file From b9e8f20ff1e4685245fce7597ce8f4406d883e01 Mon Sep 17 00:00:00 2001 From: xtalax Date: Tue, 21 Mar 2023 15:34:36 +0000 Subject: [PATCH 24/26] fix --- src/scalar_discretization.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/scalar_discretization.jl b/src/scalar_discretization.jl index 04dbbf56b..c17d1bcec 100644 --- a/src/scalar_discretization.jl +++ b/src/scalar_discretization.jl @@ -52,7 +52,6 @@ function generate_system(alleqs, bceqs, ics, discvars, defaults, ps, tspan, meta else checks = true end - end try if t === nothing # At the time of writing, NonlinearProblems require that the system of equations be in this form: From 49fd050a98aec84211425a2f75eed38c41d69ee8 Mon Sep 17 00:00:00 2001 From: xtalax Date: Tue, 21 Mar 2023 15:35:09 +0000 Subject: [PATCH 25/26] fix --- src/scalar_discretization.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/scalar_discretization.jl b/src/scalar_discretization.jl index c17d1bcec..59dba6d3d 100644 --- a/src/scalar_discretization.jl +++ b/src/scalar_discretization.jl @@ -47,11 +47,11 @@ function generate_system(alleqs, bceqs, ics, discvars, defaults, ps, tspan, meta alleqs = vcat(alleqs, unique(bceqs)) alldepvarsdisc = vec(reduce(vcat, vec(unique(reduce(vcat, vec.(values(discvars))))))) # Finalize - if haskey(metadata.disc.kwargs, :checks) - checks = metadata.disc.kwargs[:checks] - else - checks = true - end + # if haskey(metadata.disc.kwargs, :checks) + # checks = metadata.disc.kwargs[:checks] + # else + # checks = true + # end try if t === nothing # At the time of writing, NonlinearProblems require that the system of equations be in this form: From 4a496ffc05afebf29b9e4b081ef0dc2fbe100758 Mon Sep 17 00:00:00 2001 From: xtalax Date: Tue, 21 Mar 2023 15:35:41 +0000 Subject: [PATCH 26/26] fix --- src/scalar_discretization.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/scalar_discretization.jl b/src/scalar_discretization.jl index 59dba6d3d..bc5e7d207 100644 --- a/src/scalar_discretization.jl +++ b/src/scalar_discretization.jl @@ -50,7 +50,7 @@ function generate_system(alleqs, bceqs, ics, discvars, defaults, ps, tspan, meta # if haskey(metadata.disc.kwargs, :checks) # checks = metadata.disc.kwargs[:checks] # else - # checks = true + checks = true # end try if t === nothing