From f770ff15baf29c3b9a4316c5da82e16220c14246 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Thu, 20 Jul 2023 14:43:18 +0100 Subject: [PATCH 01/48] fix halfar --- src/system_parsing/pde_system_transformation.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/system_parsing/pde_system_transformation.jl b/src/system_parsing/pde_system_transformation.jl index 1328ef22c..556019b4a 100644 --- a/src/system_parsing/pde_system_transformation.jl +++ b/src/system_parsing/pde_system_transformation.jl @@ -9,7 +9,7 @@ function PDEBase.transform_pde_system!(v::PDEBase.VariableMap, boundarymap, sys: eqs = copy(sys.eqs) bcs = copy(sys.bcs) done = false - # Replace bad terms with a greedy strategy until the system comes up clean + # Replace bad terms with a greedy strategy until the system comes up clean. while !done done = true for eq in eqs @@ -103,6 +103,10 @@ function nonlinlap_check(term, differential) end end + if length(findall(has_derivatives, args)) > 1 + return nothing + end + is = findall(args) do arg if istree(arg) op = operation(arg) From dcf9b9fe1dd220b4c593fc397040dc5315c8317e Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Thu, 20 Jul 2023 16:32:21 +0100 Subject: [PATCH 02/48] Beef up nonlinlap handling --- .../nonlinear_laplacian.jl | 58 ++++++++++++++++--- .../pde_system_transformation.jl | 4 -- 2 files changed, 49 insertions(+), 13 deletions(-) diff --git a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl index 0627af48c..c7862c8ad 100644 --- a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl +++ b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl @@ -26,9 +26,10 @@ where `finitediff(u, i)` is the finite difference at the interpolated point `i` And so on. """ -function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace, bs, depvars, x, u) +function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace, bcmap, depvars, x, u) # Based on the paper https://web.mit.edu/braatzgroup/analysis_of_finite_difference_discretization_schemes_for_diffusion_in_spheres_with_variable_diffusivity.pdf # See scheme 1, namely the term without the 1/r dependence. See also #354 and #371 in DiffEqOperators, the previous home of this package. + bs = filter_interfaces(bcmap[operation(u)][x]) N = ndims(u, s) N == 0 && return 0 jx = j, x = (x2i(s, u, x), x) @@ -47,8 +48,8 @@ function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace, interface_wrap(stencil) = bwrap.(stencil, (bs,), (s,), (jx,)) # Get the correct weights and stencils for this II - inner_deriv_weights_and_stencil = [get_half_offset_weights_and_stencil(D_inner, I, s, bs, u, jx) for I in outerstencil] interp_weights_and_stencil = [get_half_offset_weights_and_stencil(inner_interpolater, I, s, bs, u, jx) for I in outerstencil] + deriv_weights_and_stencil(u, i, order) = get_half_offset_weights_and_stencil(derivweights.halfoffsetmap[1][Differential(x)^order], outerstencil[i], s, bs, u, x2i(s, u, x)) for D in (D_inner, D_outer) # map variables to symbolically inerpolated/extrapolated expressions map_vars_to_interpolated(stencil, weights) = [v => sym_dot(weights, s.discvars[v][interface_wrap(stencil)]) for v in depvars] @@ -57,35 +58,74 @@ function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace, #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])]) + # Go ham and try to discretize anything that appears inside the nonlinear laplacian + drules = generate_deriv_rules(II, s, depvars, derivweights, bmap, jx, outerstencil, deriv_weights_and_stencil) + # Take the inner finite difference 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 - interpolated_expr = map(interp_weights_and_stencil) do (weights, stencil) - substitute(substitute(expr, map_vars_to_interpolated(stencil, weights)), map_ivs_to_interpolated(stencil, weights)) + interpolated_expr = map(enumerate(interp_weights_and_stencil)) do (i, (weights, stencil)) + rules = vcat(drules(i), map_vars_to_interpolated(stencil, weights), map_ivs_to_interpolated(stencil, weights)) + substitute(expr*Differential(x)(u), rules) end # multiply the inner finite difference by the interpolated expression, and finally take the outer finite difference return sym_dot(outerweights, inner_difference .* interpolated_expr) end +function generate_deriv_rules(II, s::DiscreteSpace, depvars, derivweights, bmap, jx, outerstencil, deriv_weights_and_stencil) + # Generate rules for the derivatives of the variables + j, x = jx + mapreduce(vcat, depvars) do u + map(ivs(u, s)) do y + let orders = derivweights.orders[x] + mapreduce(vcat, orders) do order + # If we're differentiating with respect to the same variable, we need to use the correct weights and stencil + # for the order of the derivative + if isequal(x, y) + (i) -> begin + let weights, stencil = deriv_weights_and_stencil(u, i, order) + (Differential(x)^order)(u) => sym_dot(weights, s.discvars[u][interface_wrap(stencil)]) + end + end + else + # Otherwise, we will use the usual rules shifted + (i) -> begin + let II = outerstencil[i] + central_deriv_rules_cartesian = generate_cartesian_rules(II, s, depvars, derivweights, bmap, indexmap, nothing) + advection_rules = generate_advection_rules(derivweights.advection_scheme, II, s, depvars, derivweights, bmap, indexmap, nothing) + advection_rules = vcat(advection_rules, + generate_winding_rules(II, s, depvars, derivweights, bmap, + indexmap, nothing; skip = [1])) + vcat(central_deriv_rules_cartesian, advection_rules) + end + end + end + end + end + end +end + + + function replacevals(ex, s, u, depvars, II, indexmap) rules = valmaps(s, u, depvars, II, indexmap) substitute(ex, rules) 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 ivs(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, bcmap, 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 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, bcmap, 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 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, bcmap, 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 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, bcmap, 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))(*(~~a, $(Differential(x))(u), ~~d))), ~~c), ~e) => /(*(~b..., ~c..., cartesian_nonlinear_laplacian(*(~a..., ~d...), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u)), replacevals(~e, s, u, depvars, II, indexmap)) for x in ivs(u, s)]) for u in depvars], init = [])) + rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule /(*(~~b, ($(Differential(x))(*(~~a, $(Differential(x))(u), ~~d))), ~~c), ~e) => /(*(~b..., ~c..., cartesian_nonlinear_laplacian(*(~a..., ~d...), Idx(II, s, u, indexmap), derivweights, s, bcmap, depvars, x, u)), replacevals(~e, s, u, depvars, II, indexmap)) for x in ivs(u, s)]) for u in depvars], init = [])) nonlinlap_rules = [] diff --git a/src/system_parsing/pde_system_transformation.jl b/src/system_parsing/pde_system_transformation.jl index 556019b4a..4573240c8 100644 --- a/src/system_parsing/pde_system_transformation.jl +++ b/src/system_parsing/pde_system_transformation.jl @@ -103,10 +103,6 @@ function nonlinlap_check(term, differential) end end - if length(findall(has_derivatives, args)) > 1 - return nothing - end - is = findall(args) do arg if istree(arg) op = operation(arg) From e53befd8d047e7dfe65e5363b188efa11a3882fd Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Thu, 20 Jul 2023 16:35:43 +0100 Subject: [PATCH 03/48] add unit dome test, new group for adv nonlinlap --- .../nonlinear_laplacian_advanced.jl | 88 +++++++++++++++++++ test/runtests.jl | 7 ++ 2 files changed, 95 insertions(+) create mode 100644 test/pde_systems/nonlinear_laplacian_advanced.jl diff --git a/test/pde_systems/nonlinear_laplacian_advanced.jl b/test/pde_systems/nonlinear_laplacian_advanced.jl new file mode 100644 index 000000000..8b2a9300e --- /dev/null +++ b/test/pde_systems/nonlinear_laplacian_advanced.jl @@ -0,0 +1,88 @@ +using ModelingToolkit, MethodOfLines, LinearAlgebra, Test, OrdinaryDiffEq, DomainSets + +function halfar_dome(t, x, y, R0, H0, ρ, A=1e-16) + n = 3.0 + grav = 9.8101 + alpha = 1.0/9.0 + beta = 1.0/18.0 + + Gamma = 2.0/(n+2.0) * A * (ρ * grav)^n + + xcenter = 0.0 + ycenter = 0.0 + + t0 = (beta/Gamma) * (7.0/4.0)^3 * (R0^4/H0^7) # Note: this line assumes n=3! + tr=(t+t0)/t0 + + H=zeros(length(y), length(x)) + for i in eachindex(x) + for j in eachindex(y) + r = sqrt((x[i]-xcenter)^2 + (y[j]-ycenter)^2) + r=r/R0 + inside = max(0.0, 1.0 - (r / tr^beta)^((n+1.0) / n)) + + H[i,j] = H0 * inside^(n / (2.0*n+1.0)) / tr^alpha + end + end + return H + +end +@testset "Halfar ice dome glacier model." begin + + rmax = 2*1000 + rmin = -rmax + + H0 = 100 + R0 = 1000 + + asf = (t, x, y) -> halfar_dome(t, x, y, R0, H0, 917) + + @parameters x, y, t + + @variables H(..) inHx(..) inHy(..) + + Dx = Differential(x) + Dy = Differential(y) + Dt = Differential(t) + Dxx = Differential(x)^2 + Dyy = Differential(y)^2 + + n = 3.0 + grav = 9.8101 + A = 1e-16 + ρ = 910.0 + + Γ = 2.0/(n+2.0) * A * (ρ * grav)^n + + eqs = [Dt(H(t,x,y)) ~ Dx(Γ*H(t, x, y)^(n+2)*(abs(Dx(H(t, x, y)))^(n-1))*Dx(H(t, x, y))) + + Dy(Γ*H(t, x, y)^(n+2)*(abs(Dy(H(t, x, y)))^(n-1))*Dy(H(t, x, y)))] + + bcs = [H(0.0, x, y) ~ asf(0.0, x, y), + H(t, xmin, y) ~ 0.0, + H(t, xmax, y) ~ 0.0, + H(t, x, ymin) ~ 0.0, + H(t, x, ymax) ~ 0.0] + + domains = [x ∈ IntervalDomain(rmin, rmax), + y ∈ IntervalDomain(rmin, rmax), + t ∈ IntervalDomain(0.0, 1.0e6)] + + @named pdesys = PDESystem(eqs, bcs, domains, [x, y, t], [H(t, x, y)]) + + disc = MOLFiniteDifference([x => 24, y => 24], t) + + prob = discretize(pdesys, disc) + + sol = solve(prob, FBDF()) + + @test sol.retcode == :Success + + solx = sol.x + soly = sol.y + solt = sol.t + + solexact = [asf(dt, dx, dy) for dt in solt, dx in solx, dy in soly] + + @test sum(abs2, sol[H(t, x, y)] .- solexact) < 1e-2 + +end diff --git a/test/runtests.jl b/test/runtests.jl index 5916d5acf..476796b8e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -127,4 +127,11 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") include("pde_systems/MOL_1D_Linear_Convection.jl") end end + + if GROUP == "All" || GROUP == "Nonlinlap_ADV" + @time @safetestset "MOLFiniteDifference Interface: Advanced Nonlinear Diffusion" begin + include("pde_systems/nonlinear_laplacian_advanced.jl") + end + end + end end From 0d6a4b6625fb4de01dfd5bd365fcb9018365e370 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Mon, 24 Jul 2023 15:47:09 +0100 Subject: [PATCH 04/48] diffusion passing --- .../nonlinear_laplacian.jl | 61 +++++++++---------- .../spherical_laplacian.jl | 12 ++-- 2 files changed, 36 insertions(+), 37 deletions(-) diff --git a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl index c7862c8ad..e370011eb 100644 --- a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl +++ b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl @@ -49,7 +49,7 @@ function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace, # Get the correct weights and stencils for this II interp_weights_and_stencil = [get_half_offset_weights_and_stencil(inner_interpolater, I, s, bs, u, jx) for I in outerstencil] - deriv_weights_and_stencil(u, i, order) = get_half_offset_weights_and_stencil(derivweights.halfoffsetmap[1][Differential(x)^order], outerstencil[i], s, bs, u, x2i(s, u, x)) for D in (D_inner, D_outer) + deriv_weights_and_stencil(u, i, order) = get_half_offset_weights_and_stencil(derivweights.halfoffsetmap[1][Differential(x)^order], outerstencil[i], s, bs, u, (x2i(s, u, x), x)) # map variables to symbolically inerpolated/extrapolated expressions map_vars_to_interpolated(stencil, weights) = [v => sym_dot(weights, s.discvars[v][interface_wrap(stencil)]) for v in depvars] @@ -59,49 +59,46 @@ function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace, 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])]) # Go ham and try to discretize anything that appears inside the nonlinear laplacian - drules = generate_deriv_rules(II, s, depvars, derivweights, bmap, jx, outerstencil, deriv_weights_and_stencil) + drules = generate_deriv_rules(II, s, depvars, derivweights, bcmap, jx, outerstencil, deriv_weights_and_stencil, interface_wrap) - # Take the inner finite difference - 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 - - - interpolated_expr = map(enumerate(interp_weights_and_stencil)) do (i, (weights, stencil)) - rules = vcat(drules(i), map_vars_to_interpolated(stencil, weights), map_ivs_to_interpolated(stencil, weights)) + # Symbolically interpolate the expression + interpolated_expr = map(enumerate(interp_weights_and_stencil)) do (i, weights_stencil) + weights, stencil = weights_stencil + rules = vcat(mapreduce(d -> d(i), vcat, drules), map_vars_to_interpolated(stencil, weights), map_ivs_to_interpolated(stencil, weights)) substitute(expr*Differential(x)(u), rules) end # multiply the inner finite difference by the interpolated expression, and finally take the outer finite difference - return sym_dot(outerweights, inner_difference .* interpolated_expr) + return sym_dot(outerweights, interpolated_expr) end -function generate_deriv_rules(II, s::DiscreteSpace, depvars, derivweights, bmap, jx, outerstencil, deriv_weights_and_stencil) +function generate_deriv_rules(II, s::DiscreteSpace, depvars, derivweights, bmap, jx, outerstencil, deriv_weights_and_stencil, interface_wrap) # Generate rules for the derivatives of the variables j, x = jx mapreduce(vcat, depvars) do u - map(ivs(u, s)) do y + mapreduce(vcat, ivs(u, s)) do y let orders = derivweights.orders[x] - mapreduce(vcat, orders) do order - # If we're differentiating with respect to the same variable, we need to use the correct weights and stencil - # for the order of the derivative - if isequal(x, y) - (i) -> begin - let weights, stencil = deriv_weights_and_stencil(u, i, order) - (Differential(x)^order)(u) => sym_dot(weights, s.discvars[u][interface_wrap(stencil)]) + map(orders) do order + # If we're differentiating with respect to the same variable, we need to use the correct weights and stencil + # for the order of the derivative + if isequal(x, y) + (i) -> begin + let (weights, stencil) = deriv_weights_and_stencil(u, i, order) + [(Differential(x)^order)(u) => sym_dot(weights, s.discvars[u][interface_wrap(stencil)])] + end + end + else + # Otherwise, we will use the usual rules shifted + (i) -> begin + let II = outerstencil[i] + central_deriv_rules_cartesian = generate_cartesian_rules(II, s, depvars, derivweights, bmap, indexmap, nothing) + advection_rules = generate_advection_rules(derivweights.advection_scheme, II, s, depvars, derivweights, bmap, indexmap, nothing) + advection_rules = vcat(advection_rules, + generate_winding_rules(II, s, depvars, derivweights, bmap, + indexmap, nothing; skip = [1])) + vcat(central_deriv_rules_cartesian, advection_rules) + end end - end - else - # Otherwise, we will use the usual rules shifted - (i) -> begin - let II = outerstencil[i] - central_deriv_rules_cartesian = generate_cartesian_rules(II, s, depvars, derivweights, bmap, indexmap, nothing) - advection_rules = generate_advection_rules(derivweights.advection_scheme, II, s, depvars, derivweights, bmap, indexmap, nothing) - advection_rules = vcat(advection_rules, - generate_winding_rules(II, s, depvars, derivweights, bmap, - indexmap, nothing; skip = [1])) - vcat(central_deriv_rules_cartesian, advection_rules) - end end end end diff --git a/src/discretization/schemes/spherical_laplacian/spherical_laplacian.jl b/src/discretization/schemes/spherical_laplacian/spherical_laplacian.jl index c56a248df..c9bb864fb 100644 --- a/src/discretization/schemes/spherical_laplacian/spherical_laplacian.jl +++ b/src/discretization/schemes/spherical_laplacian/spherical_laplacian.jl @@ -8,8 +8,10 @@ Based on https://web.mit.edu/braatzgroup/analysis_of_finite_difference_discretiz See scheme 1 in appendix A. The r = 0 case is treated in a later appendix """ -function spherical_diffusion(innerexpr, II, derivweights, s, bs, depvars, r, u) +function spherical_diffusion(innerexpr, II, derivweights, s, bcmap, depvars, r, u) # Based on the paper https://web.mit.edu/braatzgroup/analysis_of_finite_difference_discretization_schemes_for_diffusion_in_spheres_with_variable_diffusivity.pdf + bs = filter_interfaces(bcmap[operation(u)][r]) + D_1 = derivweights.map[Differential(r)] D_2 = derivweights.map[Differential(r)^2] @@ -31,17 +33,17 @@ function spherical_diffusion(innerexpr, II, derivweights, s, bs, depvars, r, u) D_1_u = central_difference(D_1, II, s, bs, (s.x2i[r], r), u, ufunc_u) # See scheme 1 in appendix A of the paper - return exprhere*(D_1_u/substitute(r, _rsubs(r, II)) + cartesian_nonlinear_laplacian(innerexpr, II, derivweights, s, bs, depvars, r, u)) + return exprhere*(D_1_u/substitute(r, _rsubs(r, II)) + cartesian_nonlinear_laplacian(innerexpr, II, derivweights, s, bcmap, depvars, r, u)) 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...) + 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, bcmap , depvars, r, u), ~b...) 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)) + 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, bcmap, depvars, r, u)) 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) + 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, bcmap, depvars, r, u) for r in ivs(u, s)]) for u in depvars], init = [])) spherical_diffusion_rules = [] From 5eba37b147334c77d00ef13e77dac47fe8c4088e Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Mon, 31 Jul 2023 15:49:54 +0100 Subject: [PATCH 05/48] some fixes --- .../nonlinear_laplacian.jl | 24 ++++---- .../spherical_laplacian.jl | 10 ++-- src/system_parsing/interior_map.jl | 5 +- test/pde_systems/MOL_1D_Linear_Diffusion.jl | 2 +- .../nonlinear_laplacian_advanced.jl | 57 ++++++++++--------- test/runtests.jl | 27 ++++----- 6 files changed, 64 insertions(+), 61 deletions(-) diff --git a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl index e370011eb..b8e99144b 100644 --- a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl +++ b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl @@ -26,7 +26,7 @@ where `finitediff(u, i)` is the finite difference at the interpolated point `i` And so on. """ -function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace, bcmap, depvars, x, u) +function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace, indexmap, bcmap, depvars, x, u) # Based on the paper https://web.mit.edu/braatzgroup/analysis_of_finite_difference_discretization_schemes_for_diffusion_in_spheres_with_variable_diffusivity.pdf # See scheme 1, namely the term without the 1/r dependence. See also #354 and #371 in DiffEqOperators, the previous home of this package. bs = filter_interfaces(bcmap[operation(u)][x]) @@ -59,7 +59,7 @@ function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace, 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])]) # Go ham and try to discretize anything that appears inside the nonlinear laplacian - drules = generate_deriv_rules(II, s, depvars, derivweights, bcmap, jx, outerstencil, deriv_weights_and_stencil, interface_wrap) + drules = generate_deriv_rules(II, s, depvars, derivweights, indexmap, bcmap, jx, outerstencil, deriv_weights_and_stencil, interface_wrap) # Symbolically interpolate the expression interpolated_expr = map(enumerate(interp_weights_and_stencil)) do (i, weights_stencil) @@ -72,13 +72,13 @@ function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace, return sym_dot(outerweights, interpolated_expr) end -function generate_deriv_rules(II, s::DiscreteSpace, depvars, derivweights, bmap, jx, outerstencil, deriv_weights_and_stencil, interface_wrap) +function generate_deriv_rules(II, s::DiscreteSpace, depvars, derivweights, indexmap, bmap, jx, outerstencil, deriv_weights_and_stencil, interface_wrap) # Generate rules for the derivatives of the variables j, x = jx mapreduce(vcat, depvars) do u mapreduce(vcat, ivs(u, s)) do y let orders = derivweights.orders[x] - map(orders) do order + map(vcat(1, orders)) do order # If we're differentiating with respect to the same variable, we need to use the correct weights and stencil # for the order of the derivative if isequal(x, y) @@ -92,11 +92,7 @@ function generate_deriv_rules(II, s::DiscreteSpace, depvars, derivweights, bmap, (i) -> begin let II = outerstencil[i] central_deriv_rules_cartesian = generate_cartesian_rules(II, s, depvars, derivweights, bmap, indexmap, nothing) - advection_rules = generate_advection_rules(derivweights.advection_scheme, II, s, depvars, derivweights, bmap, indexmap, nothing) - advection_rules = vcat(advection_rules, - generate_winding_rules(II, s, depvars, derivweights, bmap, - indexmap, nothing; skip = [1])) - vcat(central_deriv_rules_cartesian, advection_rules) + central_deriv_rules_cartesian end end end @@ -114,15 +110,15 @@ function replacevals(ex, s, u, depvars, II, indexmap) 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, bcmap, depvars, x, u), ~d...) for x in ivs(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, indexmap, bcmap, 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, bcmap, 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))(*(~~a, $(Differential(x))(u), ~~b)) => cartesian_nonlinear_laplacian(*(~a..., ~b...), Idx(II, s, u, indexmap), derivweights, s, indexmap, bcmap, 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, bcmap, 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, indexmap, bcmap, 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, bcmap, 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, indexmapp, bcmap, 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))(*(~~a, $(Differential(x))(u), ~~d))), ~~c), ~e) => /(*(~b..., ~c..., cartesian_nonlinear_laplacian(*(~a..., ~d...), Idx(II, s, u, indexmap), derivweights, s, bcmap, depvars, x, u)), replacevals(~e, s, u, depvars, II, indexmap)) for x in ivs(u, s)]) for u in depvars], init = [])) + rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule /(*(~~b, ($(Differential(x))(*(~~a, $(Differential(x))(u), ~~d))), ~~c), ~e) => /(*(~b..., ~c..., cartesian_nonlinear_laplacian(*(~a..., ~d...), Idx(II, s, u, indexmap), derivweights, s, indexmap, bcmap, depvars, x, u)), replacevals(~e, s, u, depvars, II, indexmap)) for x in ivs(u, s)]) for u in depvars], init = [])) nonlinlap_rules = [] diff --git a/src/discretization/schemes/spherical_laplacian/spherical_laplacian.jl b/src/discretization/schemes/spherical_laplacian/spherical_laplacian.jl index c9bb864fb..840ffb0ce 100644 --- a/src/discretization/schemes/spherical_laplacian/spherical_laplacian.jl +++ b/src/discretization/schemes/spherical_laplacian/spherical_laplacian.jl @@ -8,7 +8,7 @@ Based on https://web.mit.edu/braatzgroup/analysis_of_finite_difference_discretiz See scheme 1 in appendix A. The r = 0 case is treated in a later appendix """ -function spherical_diffusion(innerexpr, II, derivweights, s, bcmap, depvars, r, u) +function spherical_diffusion(innerexpr, II, derivweights, s, indexmap, bcmap, depvars, r, u) # Based on the paper https://web.mit.edu/braatzgroup/analysis_of_finite_difference_discretization_schemes_for_diffusion_in_spheres_with_variable_diffusivity.pdf bs = filter_interfaces(bcmap[operation(u)][r]) @@ -33,17 +33,17 @@ function spherical_diffusion(innerexpr, II, derivweights, s, bcmap, depvars, r, D_1_u = central_difference(D_1, II, s, bs, (s.x2i[r], r), u, ufunc_u) # See scheme 1 in appendix A of the paper - return exprhere*(D_1_u/substitute(r, _rsubs(r, II)) + cartesian_nonlinear_laplacian(innerexpr, II, derivweights, s, bcmap, depvars, r, u)) + return exprhere*(D_1_u/substitute(r, _rsubs(r, II)) + cartesian_nonlinear_laplacian(innerexpr, II, derivweights, s, indexmap, bcmap, depvars, r, u)) 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, bcmap , depvars, r, u), ~b...) + 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, indexmap, bcmap , depvars, r, u), ~b...) 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, bcmap, depvars, r, u)) + 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, indexmap, bcmap, depvars, r, u)) 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, bcmap, depvars, r, u) + 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, indexmap, bcmap, depvars, r, u) for r in ivs(u, s)]) for u in depvars], init = [])) spherical_diffusion_rules = [] diff --git a/src/system_parsing/interior_map.jl b/src/system_parsing/interior_map.jl index bd2e2e7df..315d31874 100644 --- a/src/system_parsing/interior_map.jl +++ b/src/system_parsing/interior_map.jl @@ -95,11 +95,12 @@ function calculate_stencil_extents(s, u, discretization, orders, bcmap) # Skip if periodic in x haslower, hasupper = haslowerupper(filter_interfaces(bcmap[operation(u)][x]), x) for dorder in filter(isodd, orders[x]) + ascheme = dorder == 1 ? advection_scheme : UpwindScheme() if !haslower - lowerextents[j] = max(lowerextents[j], extent(advection_scheme, dorder)) + lowerextents[j] = max(lowerextents[j], extent(ascheme, dorder)) end if !hasupper - upperextents[j] = max(upperextents[j], extent(advection_scheme, dorder)) + upperextents[j] = max(upperextents[j], extent(ascheme, dorder)) end end end diff --git a/test/pde_systems/MOL_1D_Linear_Diffusion.jl b/test/pde_systems/MOL_1D_Linear_Diffusion.jl index 09f01484b..f4ff722b4 100644 --- a/test/pde_systems/MOL_1D_Linear_Diffusion.jl +++ b/test/pde_systems/MOL_1D_Linear_Diffusion.jl @@ -42,7 +42,7 @@ using ModelingToolkit: Differential for disc in [discretization, discretization_edge, discretization_centered, discretization_approx_order4] # Convert the PDE problem into an ODE problem - prob = discretize(pdesys, disc) + prob = discrebranch tize(pdesys, disc) # Solve ODE problem # Solve ODE problem sol = solve(prob, Tsit5(), saveat=0.1) diff --git a/test/pde_systems/nonlinear_laplacian_advanced.jl b/test/pde_systems/nonlinear_laplacian_advanced.jl index 8b2a9300e..8f90aa60c 100644 --- a/test/pde_systems/nonlinear_laplacian_advanced.jl +++ b/test/pde_systems/nonlinear_laplacian_advanced.jl @@ -1,4 +1,7 @@ using ModelingToolkit, MethodOfLines, LinearAlgebra, Test, OrdinaryDiffEq, DomainSets +using Symbolics +using Symbolics: wrap, unwrap +using SciMLBase function halfar_dome(t, x, y, R0, H0, ρ, A=1e-16) n = 3.0 @@ -14,29 +17,29 @@ function halfar_dome(t, x, y, R0, H0, ρ, A=1e-16) t0 = (beta/Gamma) * (7.0/4.0)^3 * (R0^4/H0^7) # Note: this line assumes n=3! tr=(t+t0)/t0 - H=zeros(length(y), length(x)) - for i in eachindex(x) - for j in eachindex(y) - r = sqrt((x[i]-xcenter)^2 + (y[j]-ycenter)^2) - r=r/R0 - inside = max(0.0, 1.0 - (r / tr^beta)^((n+1.0) / n)) - - H[i,j] = H0 * inside^(n / (2.0*n+1.0)) / tr^alpha - end - end - return H + r = sqrt((x-xcenter)^2 + (y-ycenter)^2) + r=r/R0 + inside = max(0.0, 1.0 - (r / tr^beta)^((n+1.0) / n)) + out = H0 * inside^(n / (2.0*n+1.0)) / tr^alpha + + return out end + +H0 = 100 +R0 = 1000 + +function asf(dt, dx, dy) + halfar_dome(dt, dx, dy, R0, H0, 917) +end + @testset "Halfar ice dome glacier model." begin rmax = 2*1000 rmin = -rmax - H0 = 100 - R0 = 1000 - - asf = (t, x, y) -> halfar_dome(t, x, y, R0, H0, 917) - + xmin = ymin = rmin + xmax = ymax = rmax @parameters x, y, t @variables H(..) inHx(..) inHy(..) @@ -54,8 +57,10 @@ end Γ = 2.0/(n+2.0) * A * (ρ * grav)^n - eqs = [Dt(H(t,x,y)) ~ Dx(Γ*H(t, x, y)^(n+2)*(abs(Dx(H(t, x, y)))^(n-1))*Dx(H(t, x, y))) - + Dy(Γ*H(t, x, y)^(n+2)*(abs(Dy(H(t, x, y)))^(n-1))*Dy(H(t, x, y)))] + abs∇H = sqrt(Dx(H(t, x, y))^2 + Dy(H(t, x, y)^2)) + + eqs = [Dt(H(t,x,y)) ~ Dx(Γ*H(t, x, y)^(n+2)*(abs∇H^(n-1))*Dx(H(t, x, y))) + + Dy(Γ*H(t, x, y)^(n+2)*(abs∇H^(n-1))*Dy(H(t, x, y)))] bcs = [H(0.0, x, y) ~ asf(0.0, x, y), H(t, xmin, y) ~ 0.0, @@ -63,9 +68,9 @@ end H(t, x, ymin) ~ 0.0, H(t, x, ymax) ~ 0.0] - domains = [x ∈ IntervalDomain(rmin, rmax), - y ∈ IntervalDomain(rmin, rmax), - t ∈ IntervalDomain(0.0, 1.0e6)] + domains = [x ∈ Interval(rmin, rmax), + y ∈ Interval(rmin, rmax), + t ∈ Interval(0.0, 1.0e6)] @named pdesys = PDESystem(eqs, bcs, domains, [x, y, t], [H(t, x, y)]) @@ -75,13 +80,13 @@ end sol = solve(prob, FBDF()) - @test sol.retcode == :Success + @test sol.retcode == SciMLBase.ReturnCode.Success - solx = sol.x - soly = sol.y - solt = sol.t + solx = sol[x] + soly = sol[y] + solt = sol[t] - solexact = [asf(dt, dx, dy) for dt in solt, dx in solx, dy in soly] + solexact = [asf(unwrap(dt), unwrap(dx), unwrap(dy)) for dt in solt, dx in solx, dy in soly] @test sum(abs2, sol[H(t, x, y)] .- solexact) < 1e-2 diff --git a/test/runtests.jl b/test/runtests.jl index 476796b8e..ff8054306 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,20 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") # Start Test Script @time begin + 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 == "Nonlinlap_ADV" + @time @safetestset "MOLFiniteDifference Interface: Advanced Nonlinear Diffusion" begin + include("pde_systems/nonlinear_laplacian_advanced.jl") + end + end + end + + if GROUP == "All" || GROUP == "Sol_Interface" @time @safetestset "MOLFiniteDifference Interface: Solution interface" begin include("components/solution_interface.jl") @@ -61,12 +75,6 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") 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 == "Higher_Order" @time @safetestset "MOLFiniteDifference Interface: 1D HigherOrder" begin include("pde_systems/MOL_1D_HigherOrder.jl") @@ -127,11 +135,4 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") include("pde_systems/MOL_1D_Linear_Convection.jl") end end - - if GROUP == "All" || GROUP == "Nonlinlap_ADV" - @time @safetestset "MOLFiniteDifference Interface: Advanced Nonlinear Diffusion" begin - include("pde_systems/nonlinear_laplacian_advanced.jl") - end - end - end end From 0c09bbb0bf93f3029c0550505dbedc94947b7319 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Thu, 3 Aug 2023 17:22:24 +0100 Subject: [PATCH 06/48] start jsc ext --- .../MethodOfLinesJuliaSimCompilerExt.jl | 52 +++++++++++++++++++ 1 file changed, 52 insertions(+) create mode 100644 ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl diff --git a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl new file mode 100644 index 000000000..0be4a8ce6 --- /dev/null +++ b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl @@ -0,0 +1,52 @@ +module MethodOfLinesJuliaSimCompilerExt + using MethodOfLines, JuliaSimCompiler, PDEBase + +function add_metadata!(m::JuliaSimCompiler.SystemInfo, meta) + m.system_type.metadata[] = meta +end + +function generate_system(disc_state::PDEBase.EquationState, s, u0, tspan, metadata, + disc::MethodOfLines.MOLFiniteDifference) + discvars = get_discvars(s) + t = get_time(disc) + name = metadata.pdesys.name + pdesys = metadata.pdesys + alleqs = vcat(disc_state.eqs, unique(disc_state.bceqs)) + alldepvarsdisc = vec(reduce(vcat, vec(unique(reduce(vcat, vec.(values(discvars))))))) + + defaults = Dict(pdesys.ps === nothing || pdesys.ps === SciMLBase.NullParameters() ? u0 : + vcat(u0, pdesys.ps)) + ps = pdesys.ps === nothing || pdesys.ps === SciMLBase.NullParameters() ? Num[] : + first.(pdesys.ps) + # Finalize + # 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: + # 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, 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 = IRSystem(alleqs, t, alldepvarsdisc, ps, defaults = defaults, + system_type = metadata, checks = checks) + return sys, tspan + end + catch e + println("The system of equations is:") + println(alleqs) + println() + println("Discretization failed, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.") + println() + rethrow(e) + end +end +end From ac9c2f4a61106aa04c4d1c086a2a4fe7bec4e803 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Fri, 4 Aug 2023 01:35:36 +0100 Subject: [PATCH 07/48] right interface --- .../MethodOfLinesJuliaSimCompilerExt.jl | 5 +++-- .../schemes/function_scheme/function_scheme.jl | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl index 0be4a8ce6..f45e6a80b 100644 --- a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl +++ b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl @@ -36,8 +36,9 @@ function generate_system(disc_state::PDEBase.EquationState, s, u0, tspan, metada 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 = IRSystem(alleqs, t, alldepvarsdisc, ps, defaults = defaults, - system_type = metadata, checks = checks) + sys = ODESystem(alleqs, t, alldepvarsdisc, ps, defaults = defaults, + checks = checks) + sys = IRSystem(sys) return sys, tspan end catch e diff --git a/src/discretization/schemes/function_scheme/function_scheme.jl b/src/discretization/schemes/function_scheme/function_scheme.jl index 757289e72..686d0095f 100644 --- a/src/discretization/schemes/function_scheme/function_scheme.jl +++ b/src/discretization/schemes/function_scheme/function_scheme.jl @@ -35,7 +35,7 @@ function function_scheme(F::FunctionalScheme, II, s, bs, jx, u, ufunc) # 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 = vcat(F.ps, params(s)) - t = s.time + t = Num(s.time) itap = map(I -> I[j], Itap) discx = @view s.grid[x][itap] dx = s.dxs[x] From c7542e98d4e6dec49f58d4ce35624e3613628324 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Fri, 4 Aug 2023 14:47:51 +0100 Subject: [PATCH 08/48] Ext loading --- Project.toml | 13 ++++++++++--- .../MethodOfLinesJuliaSimCompilerExt.jl | 11 +++++++---- test/pde_systems/MOL_Mixed_Deriv.jl | 2 +- test/pde_systems/MOLtest1.jl | 1 + test/pde_systems/MOLtest2.jl | 1 + test/pde_systems/brusselator_eq.jl | 2 +- test/runtests.jl | 12 ++++++------ 7 files changed, 27 insertions(+), 15 deletions(-) diff --git a/Project.toml b/Project.toml index 23ffbce7f..79f8003f2 100644 --- a/Project.toml +++ b/Project.toml @@ -20,6 +20,12 @@ SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" TermInterface = "8ea1fca8-c5ef-4a55-8b96-4e9afe9c9a3c" +[weakdeps] +JuliaSimCompiler = "8391cb6b-4921-5777-4e45-fd9aab8cb88d" + +[extensions] +MethodOfLinesJuliaSimCompilerExt = ["JuliaSimCompiler"] + [compat] Combinatorics = "1" DiffEqBase = "6" @@ -29,9 +35,9 @@ Interpolations = "0.14" Latexify = "0.15, 0.16" ModelingToolkit = "8" OrdinaryDiffEq = "6" -PDEBase = "0.1.4" SafeTestsets = "0.0.1" -SciMLBase = "1.90, 1.91" +SciMLBase = "1.94" +PDEBase = "0.1.4" StaticArrays = "1" SymbolicUtils = "1" Symbolics = "5" @@ -42,7 +48,8 @@ julia = "1.6" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" +JuliaSimCompiler = "8391cb6b-4921-5777-4e45-fd9aab8cb88d" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "NonlinearSolve", "SafeTestsets", "StableRNGs"] +test = ["Test", "NonlinearSolve", "SafeTestsets", "StableRNGs", "JuliaSimCompiler"] diff --git a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl index f45e6a80b..3034b6607 100644 --- a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl +++ b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl @@ -1,12 +1,13 @@ module MethodOfLinesJuliaSimCompilerExt using MethodOfLines, JuliaSimCompiler, PDEBase -function add_metadata!(m::JuliaSimCompiler.SystemInfo, meta) - m.system_type.metadata[] = meta +function add_metadata!(sys::JuliaSimCompiler.IRSystem, meta) + sys.info.parent.metadata.metadata[] = meta end function generate_system(disc_state::PDEBase.EquationState, s, u0, tspan, metadata, disc::MethodOfLines.MOLFiniteDifference) + println("JuliaSimCompiler: generate_system") discvars = get_discvars(s) t = get_time(disc) name = metadata.pdesys.name @@ -35,9 +36,9 @@ function generate_system(disc_state::PDEBase.EquationState, s, u0, tspan, metada 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, - checks = checks) + name = name, + metadata = metadata, checks = checks) sys = IRSystem(sys) return sys, tspan end @@ -50,4 +51,6 @@ function generate_system(disc_state::PDEBase.EquationState, s, u0, tspan, metada rethrow(e) end end + +export generate_system end diff --git a/test/pde_systems/MOL_Mixed_Deriv.jl b/test/pde_systems/MOL_Mixed_Deriv.jl index 82c4b5261..b309e39b1 100644 --- a/test/pde_systems/MOL_Mixed_Deriv.jl +++ b/test/pde_systems/MOL_Mixed_Deriv.jl @@ -1,4 +1,4 @@ -using ModelingToolkit, MethodOfLines, LinearAlgebra, Test, OrdinaryDiffEq, DomainSets +using ModelingToolkit, MethodOfLines, LinearAlgebra, Test, OrdinaryDiffEq, DomainSets, JuliaSimCompiler using ModelingToolkit: Differential # Broken in MTK diff --git a/test/pde_systems/MOLtest1.jl b/test/pde_systems/MOLtest1.jl index d237b56a9..6e31c3a71 100644 --- a/test/pde_systems/MOLtest1.jl +++ b/test/pde_systems/MOLtest1.jl @@ -3,6 +3,7 @@ using ModelingToolkit: operation, istree, arguments using DomainSets using NonlinearSolve using StableRNGs +using JuliaSimCompiler using Test # # Define some variables diff --git a/test/pde_systems/MOLtest2.jl b/test/pde_systems/MOLtest2.jl index 46ad22996..23c7e3f2c 100644 --- a/test/pde_systems/MOLtest2.jl +++ b/test/pde_systems/MOLtest2.jl @@ -3,6 +3,7 @@ using ModelingToolkit: operation, istree, arguments using DomainSets using NonlinearSolve using StableRNGs +using JuliaSimCompiler using Test # # Define some variables diff --git a/test/pde_systems/brusselator_eq.jl b/test/pde_systems/brusselator_eq.jl index 57eaf9bfd..3a78d310c 100644 --- a/test/pde_systems/brusselator_eq.jl +++ b/test/pde_systems/brusselator_eq.jl @@ -1,5 +1,5 @@ using ModelingToolkit, MethodOfLines, LinearAlgebra, OrdinaryDiffEq -using DomainSets +using DomainSets, JuliaSimCompiler # using Plots diff --git a/test/runtests.jl b/test/runtests.jl index 5916d5acf..53012a164 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,12 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") # Start Test Script @time begin + if GROUP == "All" || GROUP == "Brusselator" + @time @safetestset "MOLFiniteDifference Interface: 2D Brusselator Equation" begin + include("pde_systems/brusselator_eq.jl") + end + end + if GROUP == "All" || GROUP == "Sol_Interface" @time @safetestset "MOLFiniteDifference Interface: Solution interface" begin include("components/solution_interface.jl") @@ -110,12 +116,6 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") end end - if GROUP == "All" || GROUP == "Brusselator" - @time @safetestset "MOLFiniteDifference Interface: 2D Brusselator Equation" begin - include("pde_systems/brusselator_eq.jl") - end - end - if GROUP == "All" || GROUP == "Burgers" @time @safetestset "MOLFiniteDifference Interface: 2D Burger's Equation" begin include("pde_systems/burgers_eq.jl") From c2651e9e3753f62c45f3522314bd267ec9e1a826 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Fri, 4 Aug 2023 14:50:41 +0100 Subject: [PATCH 09/48] format --- .../MethodOfLinesJuliaSimCompilerExt.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl index 3034b6607..9b52cbd5b 100644 --- a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl +++ b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl @@ -36,8 +36,7 @@ function generate_system(disc_state::PDEBase.EquationState, s, u0, tspan, metada 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, + sys = ODESystem(alleqs, t, alldepvarsdisc, ps, defaults = defaults, name = name, metadata = metadata, checks = checks) sys = IRSystem(sys) return sys, tspan From 7aabd9420d4a7de09341aabdf1cb004d1004f4da Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Fri, 4 Aug 2023 18:17:10 +0100 Subject: [PATCH 10/48] fixes --- .../MethodOfLinesJuliaSimCompilerExt.jl | 7 +++---- test/pde_systems/brusselator_eq.jl | 8 ++++---- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl index 9b52cbd5b..0c71c734f 100644 --- a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl +++ b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl @@ -1,11 +1,11 @@ module MethodOfLinesJuliaSimCompilerExt - using MethodOfLines, JuliaSimCompiler, PDEBase - + using JuliaSimCompiler, MethodOfLines, PDEBase, SciMLBase, ModelingToolkit + using PDEBase: get_discvars, get_time function add_metadata!(sys::JuliaSimCompiler.IRSystem, meta) sys.info.parent.metadata.metadata[] = meta end -function generate_system(disc_state::PDEBase.EquationState, s, u0, tspan, metadata, +function PDEBase.generate_system(disc_state::PDEBase.EquationState, s, u0, tspan, metadata, disc::MethodOfLines.MOLFiniteDifference) println("JuliaSimCompiler: generate_system") discvars = get_discvars(s) @@ -51,5 +51,4 @@ function generate_system(disc_state::PDEBase.EquationState, s, u0, tspan, metada end end -export generate_system end diff --git a/test/pde_systems/brusselator_eq.jl b/test/pde_systems/brusselator_eq.jl index 3a78d310c..263be3bb7 100644 --- a/test/pde_systems/brusselator_eq.jl +++ b/test/pde_systems/brusselator_eq.jl @@ -1,5 +1,5 @@ -using ModelingToolkit, MethodOfLines, LinearAlgebra, OrdinaryDiffEq using DomainSets, JuliaSimCompiler +using ModelingToolkit, MethodOfLines, LinearAlgebra, OrdinaryDiffEq # using Plots @@ -7,7 +7,7 @@ using DomainSets, JuliaSimCompiler begin #@testset "Test 01: Brusselator equation 2D" begin @parameters x y t @variables u(..) v(..) - Dt = Differential(t) + Difft = Differential(t) Dx = Differential(x) Dy = Differential(y) Dxx = Differential(x)^2 @@ -26,8 +26,8 @@ begin #@testset "Test 01: Brusselator equation 2D" begin u0(x,y,t) = 22(y*(1-y))^(3/2) v0(x,y,t) = 27(x*(1-x))^(3/2) - eq = [Dt(u(x,y,t)) ~ 1. + v(x,y,t)*u(x,y,t)^2 - 4.4*u(x,y,t) + α*∇²(u(x,y,t)) + brusselator_f(x, y, t), - Dt(v(x,y,t)) ~ 3.4*u(x,y,t) - v(x,y,t)*u(x,y,t)^2 + α*∇²(v(x,y,t))] + eq = [Difft(u(x,y,t)) ~ 1. + v(x,y,t)*u(x,y,t)^2 - 4.4*u(x,y,t) + α*∇²(u(x,y,t)) + brusselator_f(x, y, t), + Difft(v(x,y,t)) ~ 3.4*u(x,y,t) - v(x,y,t)*u(x,y,t)^2 + α*∇²(v(x,y,t))] domains = [x ∈ Interval(x_min, x_max), y ∈ Interval(y_min, y_max), From c358ce4e29b5c115d324a615d199a11c4bd6594b Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Mon, 7 Aug 2023 18:17:10 +0100 Subject: [PATCH 11/48] fixes --- .../MethodOfLinesJuliaSimCompilerExt.jl | 127 +++++++++++------- 1 file changed, 81 insertions(+), 46 deletions(-) diff --git a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl index 0c71c734f..8af5b91af 100644 --- a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl +++ b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl @@ -1,54 +1,89 @@ module MethodOfLinesJuliaSimCompilerExt - using JuliaSimCompiler, MethodOfLines, PDEBase, SciMLBase, ModelingToolkit - using PDEBase: get_discvars, get_time -function add_metadata!(sys::JuliaSimCompiler.IRSystem, meta) - sys.info.parent.metadata.metadata[] = meta +using JuliaSimCompiler, MethodOfLines, PDEBase, SciMLBase, ModelingToolkit +using PDEBase: get_discvars, get_time, add_metadata! + +function PDEBase.add_metadata!(sys::JuliaSimCompiler.IRSystem, meta) + sys.info.parent.metadata.metadata[] = meta end function PDEBase.generate_system(disc_state::PDEBase.EquationState, s, u0, tspan, metadata, - disc::MethodOfLines.MOLFiniteDifference) - println("JuliaSimCompiler: generate_system") - discvars = get_discvars(s) - t = get_time(disc) - name = metadata.pdesys.name - pdesys = metadata.pdesys - alleqs = vcat(disc_state.eqs, unique(disc_state.bceqs)) - alldepvarsdisc = vec(reduce(vcat, vec(unique(reduce(vcat, vec.(values(discvars))))))) + disc::MethodOfLines.MOLFiniteDifference) + println("JuliaSimCompiler: generate_system") + discvars = get_discvars(s) + t = get_time(disc) + name = metadata.pdesys.name + pdesys = metadata.pdesys + alleqs = vcat(disc_state.eqs, unique(disc_state.bceqs)) + alldepvarsdisc = vec(reduce(vcat, vec(unique(reduce(vcat, vec.(values(discvars))))))) + + defaults = Dict(pdesys.ps === nothing || pdesys.ps === SciMLBase.NullParameters() ? u0 : + vcat(u0, pdesys.ps)) + ps = pdesys.ps === nothing || pdesys.ps === SciMLBase.NullParameters() ? Num[] : + first.(pdesys.ps) + # Finalize + # 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: + # 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, 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, checks = checks) + @show sys.eqs + sys = IRSystem(sys) + return sys, tspan + end + catch e + println("The system of equations is:") + println(alleqs) + println() + println("Discretization failed, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.") + println() + rethrow(e) + end +end + +function SciMLBase.discretize(pdesys::PDESystem, + discretization::MOLFiniteDifference; + analytic = nothing, kwargs...) + sys, tspan = SciMLBase.symbolic_discretize(pdesys, discretization) + try + simpsys = structural_simplify(sys) + if tspan === nothing + return prob = NonlinearProblem(simpsys, ones(length(simpsys.states)); + discretization.kwargs..., kwargs...) + else + # Use ODAE if nessesary + + prob = ODEProblem(simpsys, Pair[], tspan; discretization.kwargs..., + kwargs...) + if analytic === nothing + return prob + else + f = ODEFunction(pdesys, discretization, analytic = analytic, + discretization.kwargs..., kwargs...) - defaults = Dict(pdesys.ps === nothing || pdesys.ps === SciMLBase.NullParameters() ? u0 : - vcat(u0, pdesys.ps)) - ps = pdesys.ps === nothing || pdesys.ps === SciMLBase.NullParameters() ? Num[] : - first.(pdesys.ps) - # Finalize - # 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: - # 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, 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, checks = checks) - sys = IRSystem(sys) - return sys, tspan - end - catch e - println("The system of equations is:") - println(alleqs) - println() - println("Discretization failed, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.") - println() - rethrow(e) - end + return ODEProblem(f, prob.u0, prob.tspan, prob.p; + discretization.kwargs..., kwargs...) + end + end + catch e + PDEBase.error_analysis(sys, e) + end end +function PDEBase.error_analysis(sys::IRSystem, e) + rethrow(e) + end +end \ No newline at end of file From 08f74f419029b972b515727c6e7b62f7e889114b Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Wed, 30 Aug 2023 18:24:18 +0100 Subject: [PATCH 12/48] no show --- .../MethodOfLinesJuliaSimCompilerExt.jl | 1 + test/pde_systems/MOL_1D_Linear_Diffusion.jl | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl index 8af5b91af..b1565daac 100644 --- a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl +++ b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl @@ -40,6 +40,7 @@ function PDEBase.generate_system(disc_state::PDEBase.EquationState, s, u0, tspan sys = ODESystem(alleqs, t, alldepvarsdisc, ps, defaults = defaults, name = name, metadata = metadata, checks = checks) @show sys.eqs + @show sys.observed sys = IRSystem(sys) return sys, tspan end diff --git a/test/pde_systems/MOL_1D_Linear_Diffusion.jl b/test/pde_systems/MOL_1D_Linear_Diffusion.jl index 09f01484b..67b599c55 100644 --- a/test/pde_systems/MOL_1D_Linear_Diffusion.jl +++ b/test/pde_systems/MOL_1D_Linear_Diffusion.jl @@ -1,7 +1,7 @@ # 1D diffusion problem # Packages and inclusions -using ModelingToolkit, MethodOfLines, LinearAlgebra, Test, OrdinaryDiffEq, DomainSets +using ModelingToolkit, MethodOfLines, LinearAlgebra, Test, OrdinaryDiffEq, DomainSets, JuliaSimCompiler using ModelingToolkit: Differential # Tests From a59a44bb76d74a1c3c617cd7a00369d251255927 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Mon, 4 Sep 2023 13:58:15 +0100 Subject: [PATCH 13/48] update tests --- .../MethodOfLinesJuliaSimCompilerExt.jl | 2 - test/pde_systems/MOL_1D_Linear_Diffusion.jl | 2 +- test/pde_systems/MOL_Mixed_Deriv.jl | 2 +- test/pde_systems/MOLtest1.jl | 1 - test/pde_systems/MOLtest1_JSC.jl | 353 +++++++++++++ test/pde_systems/MOLtest2.jl | 1 - test/pde_systems/MOLtest2_JSC.jl | 478 ++++++++++++++++++ test/runtests.jl | 2 + 8 files changed, 835 insertions(+), 6 deletions(-) create mode 100644 test/pde_systems/MOLtest1_JSC.jl create mode 100644 test/pde_systems/MOLtest2_JSC.jl diff --git a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl index b1565daac..0d58c23c0 100644 --- a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl +++ b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl @@ -39,8 +39,6 @@ function PDEBase.generate_system(disc_state::PDEBase.EquationState, s, u0, tspan # * 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, checks = checks) - @show sys.eqs - @show sys.observed sys = IRSystem(sys) return sys, tspan end diff --git a/test/pde_systems/MOL_1D_Linear_Diffusion.jl b/test/pde_systems/MOL_1D_Linear_Diffusion.jl index 67b599c55..09f01484b 100644 --- a/test/pde_systems/MOL_1D_Linear_Diffusion.jl +++ b/test/pde_systems/MOL_1D_Linear_Diffusion.jl @@ -1,7 +1,7 @@ # 1D diffusion problem # Packages and inclusions -using ModelingToolkit, MethodOfLines, LinearAlgebra, Test, OrdinaryDiffEq, DomainSets, JuliaSimCompiler +using ModelingToolkit, MethodOfLines, LinearAlgebra, Test, OrdinaryDiffEq, DomainSets using ModelingToolkit: Differential # Tests diff --git a/test/pde_systems/MOL_Mixed_Deriv.jl b/test/pde_systems/MOL_Mixed_Deriv.jl index b309e39b1..82c4b5261 100644 --- a/test/pde_systems/MOL_Mixed_Deriv.jl +++ b/test/pde_systems/MOL_Mixed_Deriv.jl @@ -1,4 +1,4 @@ -using ModelingToolkit, MethodOfLines, LinearAlgebra, Test, OrdinaryDiffEq, DomainSets, JuliaSimCompiler +using ModelingToolkit, MethodOfLines, LinearAlgebra, Test, OrdinaryDiffEq, DomainSets using ModelingToolkit: Differential # Broken in MTK diff --git a/test/pde_systems/MOLtest1.jl b/test/pde_systems/MOLtest1.jl index 6e31c3a71..d237b56a9 100644 --- a/test/pde_systems/MOLtest1.jl +++ b/test/pde_systems/MOLtest1.jl @@ -3,7 +3,6 @@ using ModelingToolkit: operation, istree, arguments using DomainSets using NonlinearSolve using StableRNGs -using JuliaSimCompiler using Test # # Define some variables diff --git a/test/pde_systems/MOLtest1_JSC.jl b/test/pde_systems/MOLtest1_JSC.jl new file mode 100644 index 000000000..6e31c3a71 --- /dev/null +++ b/test/pde_systems/MOLtest1_JSC.jl @@ -0,0 +1,353 @@ +using ModelingToolkit, MethodOfLines, LinearAlgebra, OrdinaryDiffEq +using ModelingToolkit: operation, istree, arguments +using DomainSets +using NonlinearSolve +using StableRNGs +using JuliaSimCompiler +using Test + +# # Define some variables +@testset "Heat Equation 1D 2 variables" begin + @parameters t x + @variables u(..) v(..) + Dt = Differential(t) + Dx = Differential(x) + Dxx = Differential(x)^2 + eqs = [Dt(u(t, x)) ~ Dxx(u(t, x)), + Dt(v(t, x)) ~ Dxx(v(t, x))] + bcs = [u(0, x) ~ -x * (x - 1) * sin(x), + v(0, x) ~ -x * (x - 1) * sin(x), + u(t, 0) ~ 0.0, u(t, 1) ~ 0.0, + v(t, 0) ~ 0.0, v(t, 1) ~ 0.0] + + domains = [t ∈ Interval(0.0, 1.0), + x ∈ Interval(0.0, 1.0)] + + @named pdesys = PDESystem(eqs, bcs, domains, [t, x], [u(t, x), v(t, x)]) + discretization = MOLFiniteDifference([x => 0.1], t; grid_align=edge_align) + prob = discretize(pdesys, discretization) # This gives an ODEProblem since it's time-dependent + sol = solve(prob, Tsit5()) +end + +@testset "Heat Equation 2D 1 variable" begin + @parameters t x y + @variables u(..) + Dxx = Differential(x)^2 + Dyy = Differential(y)^2 + Dt = Differential(t) + t_min = 0.0 + t_max = 2.0 + x_min = 0.0 + x_max = 2.0 + y_min = 0.0 + y_max = 2.0 + + # 3D PDE + eq = Dt(u(t, x, y)) ~ Dxx(u(t, x, y)) + Dyy(u(t, x, y)) + + analytic_sol_func(t, x, y) = exp(x + y) * cos(x + y + 4t) + # Initial and boundary conditions + bcs = [u(t_min, x, y) ~ analytic_sol_func(t_min, x, y), + u(t, x_min, y) ~ analytic_sol_func(t, x_min, y), + u(t, x_max, y) ~ analytic_sol_func(t, x_max, y), + u(t, x, y_min) ~ analytic_sol_func(t, x, y_min), + u(t, x, y_max) ~ analytic_sol_func(t, x, y_max)] + + # Space and time domains + domains = [t ∈ Interval(t_min, t_max), + x ∈ Interval(x_min, x_max), + y ∈ Interval(y_min, y_max)] + @named pdesys = PDESystem([eq], bcs, domains, [t, x, y], [u(t, x, y)]) + # Method of lines discretization + dx = 0.1 + dy = 0.2 + discretization = MOLFiniteDifference([x => dx, y => dy], t) + prob = ModelingToolkit.discretize(pdesys, discretization) + sol = solve(prob, Tsit5()) +end + +# Diffusion in a sphere +@testset "Spherical Diffusion" begin + + @parameters t r + @variables u(..) + Dt = Differential(t) + Dr = Differential(r) + Drr = Dr^2 + eq = Dt(u(t, r)) ~ 1 / r^2 * Dr(r^2 * Dr(u(t, r))) + bcs = [u(0, r) ~ -r * (r - 1) * sin(r), + Dr(u(t, 0)) ~ 0.0, u(t, 1) ~ sin(1)] + + domains = [t ∈ Interval(0.0, 1.0), + r ∈ Interval(0.0, 1.0)] + + @named pdesys = PDESystem(eq, bcs, domains, [t, r], [u(t, r)]) + discretization = MOLFiniteDifference([r => 0.1], t) + prob = discretize(pdesys, discretization) # This gives an ODEProblem since it's time-dependent + sol = solve(prob, Tsit5()) +end + +@testset "RHS = 0" begin + @parameters x, t + @variables u(..) + Dt = Differential(t) + Dx = Differential(x) + Dx2 = Differential(x)^2 + Dx3 = Differential(x)^3 + Dx4 = Differential(x)^4 + + α = 1 + β = 4 + γ = 1 + eq = Dt(u(x, t)) + u(x, t) * Dx(u(x, t)) + α * Dx2(u(x, t)) + β * Dx3(u(x, t)) + γ * Dx4(u(x, t)) ~ 0 + + du(x, t; z=-x / 2 + t) = 15 / 2 * (tanh(z) + 1) * (3 * tanh(z) - 1) * sech(z)^2 + + bcs = [u(x, 0) ~ x^2, + Dx(u(-10, t)) ~ du(-10, t), + Dx(u(10, t)) ~ du(10, t)] + + # Space and time domains + domains = [x ∈ Interval(-10.0, 10.0), + t ∈ Interval(0.0, 0.5)] + # Discretization + dx = 0.4 + dt = 0.2 + + discretization = MOLFiniteDifference([x => dx], t; approx_order=4) + @named pdesys = PDESystem(eq, bcs, domains, [x, t], [u(x, t)]) + prob = discretize(pdesys, discretization) +end + +@testset "Wave Equation" begin + # Parameters, variables, and derivatives + @parameters t x + @variables u(..) + Dt = Differential(t) # Required in ICs and equation + Dtt = Differential(t)^2 # required in equation + Dxxxx = Differential(x)^4 # required in equation + Dxx = Differential(x)^2 # required in BCs + # some parameters + EI = 291.6667 + m = 1.3850 + c = 1.0 + p = 0.01 + L = 2.0 + + # 1D PDE and boundaru conditions + eq = m * Dtt(u(t, x)) + c * Dt(u(t, x)) + EI * Dxxxx(u(t, x)) ~ 0 + + ic_bc = [u(0, x) ~ (p * x * (x^3 + L^3 - 2 * L * x^2) / (24 * EI)), #for all 0 < u < L + Dt(u(0, x)) ~ 0.0, # for all 0 < u < L + u(t, 0) ~ 0.0, # for all t > 0,, displacement zero at u=0 + u(t, 2) ~ 0.0, # for all t > 0,, displacement zero at u=L + Dxx(u(t, 0)) ~ 0.0, # for all t > 0,, curvature zero at u=0 + Dxx(u(t, 2)) ~ 0.0] # for all t > 0,, curvature zero at u=L + + # Space and time domains + domains = [t ∈ Interval(0.0, 2.0), + x ∈ Interval(0.0, L)] + + dt = 0.1 # dt related to saving the data.. not actual dt + # PDE sustem + @named pdesys = PDESystem(eq, ic_bc, domains, [t, x], [u(t, x)]) + + # Method of lines discretization + dx = 0.1 + order = 2 + discretization = MOLFiniteDifference([x => dx, t => dt], approx_order=order) + + # Convert the PDE problem into an ODE problem + prob = discretize(pdesys, discretization) + + # Solve the ODE problem + sol = NonlinearSolve.solve(prob, NewtonRaphson()) +end + +@testset "Rearranged Robin" begin + @parameters x t + + @variables c(..) + + ∂t = Differential(t) + ∂x = Differential(x) + ∂²x = Differential(x)^2 + + D₀ = 1.0 + R = 0.5 + cₑ = 2.5 + ℓ = 2.0 + α = 1 / 3 + + diff_eq = ∂t(c(x, t)) ~ ∂x((D₀ + α * c(x, t)) * ∂x(c(x, t))) + + bcs = [c(x, 0) ~ 0.0, # initial condition + R * (D₀ + α * c(0, t)) * ∂x(c(0, t)) ~ c(0, t) - cₑ, # Robin + ∂x(c(ℓ, t)) ~ 0.0] # no flux + + domains = [t ∈ Interval(0.0, 10.0), + x ∈ Interval(0.0, ℓ)] + + @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) + + # Method of lines discretization + dx = 0.1 + order = 2 + discretization = MOLFiniteDifference([x => dx], t, approx_order=order) + + # Convert the PDE problem into an ODE problem + prob = discretize(pdesys, discretization) + + sol = solve(prob, Rodas4()) +end + +@testset "Array u" begin + # Dependencies + N = 6 # number of dependent variables + + # Variables, parameters, and derivatives + @parameters x + @variables u[1:N](..) + Dx = Differential(x) + Dxx = Differential(x)^2 + + # Domain edges + x_min = 0.0 + x_max = 1.0 + + # Discretization parameters + dx = 0.1 + order = 2 + + #u = collect(u) + + # Equations + eqs = Vector{ModelingToolkit.Equation}(undef, N) + for i = 1:N + eqs[i] = Dxx(u[i](x)) ~ u[i](x) + end + + # Initial and boundary conditions + bcs = Vector{ModelingToolkit.Equation}(undef, 2 * N) + for i = 1:N + bcs[i] = Dx(u[i](x_min)) ~ 0.0 + end + + for i = 1:N + bcs[i+N] = u[i](x_max) ~ rand(StableRNG(0),) + end + + # Space and time domains + domains = [x ∈ Interval(x_min, x_max)] + + # PDE system + @named pdesys = PDESystem(eqs, bcs, domains, [x], collect([u[i](x) for i = 1:N])) + + # Method of lines discretization + discretization = MOLFiniteDifference([x => dx], nothing, approx_order=order) + prob = ModelingToolkit.discretize(pdesys, discretization) + + # # Solution of the ODE system + sol = NonlinearSolve.solve(prob, NewtonRaphson()) +end + +@testset "Non-Uniform Grid" begin + @parameters x y t + @variables u(..) v(..) + Dt = Differential(t) + Dx = Differential(x) + Dy = Differential(y) + Dxx = Differential(x)^2 + Dyy = Differential(y)^2 + + ∇²(u) = Dxx(u) + Dyy(u) + + brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.0 + + x_min = y_min = t_min = 0.0 + x_max = y_max = 1.0 + t_max = 11.5 + + α = 10.0 + + u0(x, y, t) = 22(y * (1 - y))^(3 / 2) + v0(x, y, t) = 27(x * (1 - x))^(3 / 2) + + eq = [Dt(u(x, y, t)) ~ 1.0 + Dy(v(x, y, t)) * u(x, y, t)^2 - 4.4 * u(x, y, t) + α * ∇²(u(x, y, t)) + brusselator_f(x, y, t), + Dt(v(x, y, t)) ~ -3.4 * Dx(u(x, y, t)) - v(x, y, t) * u(x, y, t)^2 + α * ∇²(v(x, y, t))] + + domains = [x ∈ Interval(x_min, x_max), + y ∈ Interval(y_min, y_max), + t ∈ Interval(t_min, t_max)] + + bcs = [u(x, y, 0) ~ u0(x, y, 0), + u(0, y, t) ~ 0.0, + u(1, y, t) ~ 0.0, + u(x, 0, t) ~ 0.0, + u(x, 1, t) ~ 0.0, + v(x, y, 0) ~ v0(x, y, 0), + v(0, y, t) ~ 0.0, + v(1, y, t) ~ 0.0, + v(x, 0, t) ~ 0.0, + v(x, 1, t) ~ 0.0] + + @named pdesys = PDESystem(eq, bcs, domains, [x, y, t], [u(x, y, t), v(x, y, t)]) + + # Method of lines discretization + N = 5 + + dx = [0.0, cumsum(rand(StableRNG(0), 0.05:0.01:0.2, 5))...] + if dx[end] < 1.0 + push!(dx, 1.0) + end + + dy = [0.0, cumsum(rand(StableRNG(0), 0.05:0.01:0.2, 5))...] + if dy[end] < 1.0 + push!(dy, 1.0) + end + + order = 2 + + discretization = MOLFiniteDifference([x => dx, y => dy], t, approx_order=order) + + prob = discretize(pdesys, discretization) + + sol = solve(prob, TRBDF2()) + +end + +@testset "2D variable connected to 1D variable at boundary #33" begin + @parameters t x r + @variables u(..) v(..) + Dt = Differential(t) + Dx = Differential(x) + Dxx = Differential(x)^2 + Dr = Differential(r) + Drr = Differential(r)^2 + + s = u(t, x) + v(t, x, 1) + + eqs = [Dt(u(t, x)) ~ Dxx(u(t, x)) + s, + Dt(v(t, x, r)) ~ Drr(v(t, x, r))] + bcs = [u(0, x) ~ 1, + v(0, x, r) ~ 1, + Dx(u(t, 0)) ~ 0.0, Dx(u(t, 1)) ~ 0.0, + Dr(v(t, x, 0)) ~ 0.0, Dr(v(t, x, 1)) ~ s] + + domains = [t ∈ Interval(0.0, 1.0), + x ∈ Interval(0.0, 1.0), + r ∈ Interval(0.0, 1.0)] + + @named pdesys = PDESystem(eqs, bcs, domains, [t, x, r], [u(t, x), v(t, x, r)]) + + # Method of lines discretization + dx = 0.1 + dr = 0.1 + order = 2 + discretization = MOLFiniteDifference([x => dx, r => dr], t, approx_order=order) + + # Convert the PDE problem into an ODE problem + prob = discretize(pdesys, discretization) + + sol = solve(prob, Tsit5()) +end diff --git a/test/pde_systems/MOLtest2.jl b/test/pde_systems/MOLtest2.jl index 23c7e3f2c..46ad22996 100644 --- a/test/pde_systems/MOLtest2.jl +++ b/test/pde_systems/MOLtest2.jl @@ -3,7 +3,6 @@ using ModelingToolkit: operation, istree, arguments using DomainSets using NonlinearSolve using StableRNGs -using JuliaSimCompiler using Test # # Define some variables diff --git a/test/pde_systems/MOLtest2_JSC.jl b/test/pde_systems/MOLtest2_JSC.jl new file mode 100644 index 000000000..23c7e3f2c --- /dev/null +++ b/test/pde_systems/MOLtest2_JSC.jl @@ -0,0 +1,478 @@ +using ModelingToolkit, MethodOfLines, LinearAlgebra, OrdinaryDiffEq +using ModelingToolkit: operation, istree, arguments +using DomainSets +using NonlinearSolve +using StableRNGs +using JuliaSimCompiler +using Test + +# # Define some variables + + + +@testset "Testing discretization of varied systems" begin + @parameters x t + + @variables c(..) + + ∂t = Differential(t) + ∂x = Differential(x) + ∂²x = Differential(x)^2 + + D₀ = 1.5 + α = 0.15 + χ = 1.2 + R = 0.1 + cₑ = 2.0 + ℓ = 1.0 + Δx = 0.1 + + bcs = [ + # initial condition + c(x, 0) ~ 0.0, + # Robin BC + ∂x(c(0.0, t)) / (1 + exp(α * (c(0.0, t) - χ))) * R * D₀ + cₑ - c(0.0, t) ~ 0.0, + # no flux BC + ∂x(c(ℓ, t)) ~ 0.0] + + + # define space-time plane + domains = [x ∈ Interval(0.0, ℓ), t ∈ Interval(0.0, 5.0)] + + @test_broken begin #@testset "Test 01: ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))" begin + D = D₀ / (1.0 + exp(α * (c(x, t) - χ))) + diff_eq = ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t))) + @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) + discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) + end + + @testset "Test 02: ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))" begin + D = 1.0 / (1.0 + exp(α * (c(x, t) - χ))) + diff_eq = ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t))) + @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) + discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) + end + + @testset "Test 03: ∂t(c(x, t)) ~ ∂x(1.0 / (1.0/D₀ + exp(α * (c(x, t) - χ))/D₀) * ∂x(c(x, t)))" begin + diff_eq = ∂t(c(x, t)) ~ ∂x(1.0 / (1.0 / D₀ + exp(α * (c(x, t) - χ)) / D₀) * ∂x(c(x, t))) + @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) + discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) + end + + @test_broken begin #@testset "Test 04: ∂t(c(x, t)) ~ ∂x(D₀ / (1.0 + exp(α * (c(x, t) - χ))) * ∂x(c(x, t)))" begin + diff_eq = ∂t(c(x, t)) ~ ∂x(D₀ / (1.0 + exp(α * (c(x, t) - χ))) * ∂x(c(x, t))) + @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) + discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) + end + + @testset "Test 05: ∂t(c(x, t)) ~ ∂x(1/x * ∂x(c(x, t)))" begin + diff_eq = ∂t(c(x, t)) ~ ∂x(1 / x * ∂x(c(x, t))) + @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) + discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) + end + + @test_broken begin #@testset "Test 06: ∂t(c(x, t)) ~ ∂x(x*∂x(c(x, t)))/c(x,t)" begin + diff_eq = ∂t(c(x, t)) ~ ∂x(x * ∂x(c(x, t))) / c(x, t) + @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) + discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) + end + + @testset "Test 07: ∂t(c(x, t)) ~ ∂x(1/(1+c(x,t)) ∂x(c(x, t)))" begin + diff_eq = ∂t(c(x, t)) ~ ∂x(1 / (1 + c(x, t)) * ∂x(c(x, t))) + @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) + discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) + end + + @test_broken begin #@testset "Test 08: ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x,t) * ∂x(c(x, t)))" begin + diff_eq = ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x, t) * ∂x(c(x, t))) + @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) + discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) + end + + @test_broken begin #@testset "Test 09: ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x,t) * ∂x(c(x, t)))/(1+c(x,t))" begin + diff_eq = c(x, t) * ∂x(c(x, t) * ∂x(c(x, t))) / (1 + c(x, t)) ~ 0 + @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) + discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) + end + + @test_broken begin #@testset "Test 10: ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x,t) * ∂x(c(x, t)))/(1+c(x,t))" begin + diff_eq = c(x, t) * ∂x(c(x, t)^(-1) * ∂x(c(x, t))) ~ 0 + @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) + discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) + end + + @testset "Test 11: ∂t(c(x, t)) ~ ∂x(1/(1+c(x,t)^2) ∂x(c(x, t)))" begin + diff_eq = ∂t(c(x, t)) ~ ∂x(1 / (1 + c(x, t)^2) * ∂x(c(x, t))) + @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) + discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) + end + +end + +@testset "Nonlinlap with flux interface boundary condition" begin + @parameters t x1 x2 + @variables c1(..) + @variables c2(..) + Dt = Differential(t) + + Dx1 = Differential(x1) + Dxx1 = Dx1^2 + + Dx2 = Differential(x2) + Dxx2 = Dx2^2 + + D1(c) = 1 + c / 10 + D2(c) = 1 / 10 + c / 10 + + eqs = [Dt(c1(t, x1)) ~ Dx1(D1(c1(t, x1)) * Dx1(c1(t, x1))), + Dt(c2(t, x2)) ~ Dx2(D2(c2(t, x2)) * Dx2(c2(t, x2)))] + + bcs = [c1(0, x1) ~ 1 + cospi(2 * x1), + c2(0, x2) ~ 1 + cospi(2 * x2), + Dx1(c1(t, 0)) ~ 0, + c1(t, 0.5) ~ c2(t, 0.5), + -D1(c1(t, 0.5)) * Dx1(c1(t, 0.5)) ~ -D2(c2(t, 0.5)) * Dx2(c2(t, 0.5)), + Dx2(c2(t, 1)) ~ 0] + + domains = [t ∈ Interval(0.0, 0.15), + x1 ∈ Interval(0.0, 0.5), + x2 ∈ Interval(0.5, 1.0)] + + @named pdesys = PDESystem(eqs, bcs, domains, + [t, x1, x2], [c1(t, x1), c2(t, x2)]) + + l = 40 + + disc = MOLFiniteDifference([x1 => l, x2 => l], t) + + prob = discretize(pdesys, disc) + + sol = solve(prob, FBDF(), saveat=0.01) + + x1_sol = sol[x1] + x2_sol = sol[x2] + t_sol = sol[t] + solc1 = sol[c1(t, x1)] + solc2 = sol[c2(t, x2)] + + solc = hcat(solc1[:, :], solc2[:, 2:end]) + + @test sol.retcode == SciMLBase.ReturnCode.Success +end + +@testset "Another boundaries appearing in equations case" begin + + g = 9.81 + + @parameters x z t + @variables φ(..) φ̃(..) η(..) + + Dt = Differential(t) + Dx = Differential(x) + Dz = Differential(z) + Dxx = Differential(x)^2 + Dzz = Differential(z)^2 + + eqs = [Dxx(φ(t, x, z)) + Dzz(φ(t, x, z)) ~ 0, + Dt(φ̃(t, x)) ~ -g * η(t, x), + Dt(η(t, x)) ~ Dz(φ(t, x, 1.0)) + ] + + bcs = [ + φ(0, x, z) ~ 0, + φ̃(0.0, x) ~ 0.0, + η(0.0, x) ~ cos(2 * π * x), + φ(t, x, 1.0) ~ φ̃(t, x), + Dx(φ(t, 0.0, z)) ~ 0.0, + Dx(φ(t, 1.0, z)) ~ 0.0, + Dz(φ(t, x, 0.0)) ~ 0.0, + Dx(φ̃(t, 0.0)) ~ 0.0, + Dx(φ̃(t, 1.0)) ~ 0.0, + Dx(η(t, 0.0)) ~ 0.0, + Dx(η(t, 1.0)) ~ 0.0, + ] + + domains = [x ∈ Interval(0.0, 1.0), + z ∈ Interval(0.0, 1.0), + t ∈ Interval(0.0, 10.0)] + + @named pdesys = PDESystem(eqs, bcs, domains, [t, x, z], + [φ(t, x, z), φ̃(t, x), η(t, x)]) + + + dx = 0.1 + dz = 0.1 + order = 2 + + discretization = MOLFiniteDifference([x => dx, z => dz], t, + approx_order=order, + grid_align=center_align) + + println("Discretization:") + prob = discretize(pdesys, discretization) +end + +@testset "Integrals in BCs" begin + β = 0.0005 + γ = 0.25 + amin = 0.0 + amax = 40.0 + + @parameters t a + @variables S(..) I(..) R(..) + Dt = Differential(t) + Da = Differential(a) + Ia = Integral(a in DomainSets.ClosedInterval(amin, amax)) + + + eqs = [Dt(S(t)) ~ -β * S(t) * Ia(I(a, t)), + Dt(I(a, t)) + Da(I(a, t)) ~ -γ * I(a, t), + Dt(R(t)) ~ γ * Ia(I(a, t))] + + bcs = [ + S(0) ~ 990.0, + I(0, t) ~ β * S(t) * Ia(I(a, t)), + I(a, 0) ~ 10.0 / 40.0, + R(0) ~ 0.0 + ] + + domains = [t ∈ (0.0, 40.0), a ∈ (0.0, 40.0)] + + @named pde_system = PDESystem(eqs, bcs, domains, [a, t], [S(t), I(a, t), R(t)]) + + da = 40 + discretization = MOLFiniteDifference([a => da], t) + + prob = MethodOfLines.discretize(pde_system, discretization) + + sol = solve(prob, FBDF()) + +end + +@testset "Dt in BCs" begin + # Parameters, variables, and derivatives + @parameters t x + @variables u(..) + Dt = Differential(t) + Dxx = Differential(x)^2 + + # 1D PDE and boundary conditions + eq = Dt(u(t, x)) ~ Dxx(u(t, x)) + bcs = [u(0, x) ~ 20, + Dt(u(t, 0)) ~ 100, # Heat source + Dt(u(t, 1)) ~ 0] # Zero flux + + # Space and time domains + domains = [t ∈ Interval(0.0, 1.0), + x ∈ Interval(0.0, 1.0)] + + # PDE system + @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)]) + + # Method of lines discretization + dx = 0.1 + order = 2 + discretization = MOLFiniteDifference([x => dx], t) + + # Convert the PDE problem into an ODE problem + prob = discretize(pdesys, discretization) + + # Solve ODE problem + sol = solve(prob, Rodas4(), saveat=0.2) + + discrete_x = sol[x] + discrete_t = sol[t] + solu = sol[u(t, x)] # Temperature should increase with time +end + +@testset "ODE connected to PDE at boundary" begin + @variables u(..) v(..) w(..) + @parameters t, r + Dt = Differential(t) + Dr = Differential(r) + Drr = Differential(r)^2 + + R = 1.0 + k₁ = 0.1 + k₂ = 0.1 + α = 1.0 + + u0 = 0.3 + v0 = 0.1 + w0 = 0.2 + + eq = [Dt(u(r, t)) ~ α * Drr(u(r, t)), + Dt(v(t)) ~ -k₁ * u(R, t) * v(t) + k₂ * w(t), + Dt(w(t)) ~ k₁ * u(R, t) * v(t) - k₂ * w(t) + ] + + bcs = [Dr(u(0, t)) ~ 0.0, + Dr(u(R, t)) ~ (-k₁ * u(R, t) * v(t) + k₂ * w(t)) / α, + u(r, 0) ~ u0, + v(0) ~ v0, + w(0) ~ w0 + ] + + domains = [t ∈ Interval(0.0, 10.0), + r ∈ Interval(0.0, R)] + + @named pdesys = PDESystem(eq, bcs, domains, [r, t], [u(r, t), v(t), w(t)]) + + dr = 0.1 + + disc = MOLFiniteDifference([r => dr], t) + + prob = discretize(pdesys, disc) + + sol = solve(prob, Rodas4P()) + + discrete_r = sol[r] + discrete_t = sol[t] + solu = sol[u(r, t)] + solv = sol[v(t)] + solw = sol[w(t)] +end + +@testset "ODE connected to PDE" begin + + @parameters t x + @variables u(..) v(..) + Dt = Differential(t) + Dx = Differential(x) + Dxx = Dx^2 + + # 1D PDE and boundary conditions + eqs = [Dt(u(t, x)) ~ Dxx(u(t, x)) + v(t), # This is the only line that is significantly changed from the test. + Dt(v(t)) ~ -v(t)] + + bcs = [u(0, x) ~ sin(x), + v(0) ~ 1, + u(t, 0) ~ 0, + Dx(u(t, 1)) ~ exp(-t) * cos(1)] + domains = [t ∈ Interval(0.0, 1.0), + x ∈ Interval(0.0, 1.0)] + @named pdesys = PDESystem(eqs, bcs, domains, [t, x], [u(t, x), v(t)]) + discretization = MOLFiniteDifference([x => 0.01], t) + prob = discretize(pdesys, discretization) + + sol = solve(prob, Tsit5()) + + discrete_x = sol[x] + discrete_t = sol[t] + solu = sol[u(t, x)] + solv = sol[v(t)] +end + +@testset "New style array variable conversion and interception" begin + # Parameters, variables, and derivatives + n_comp = 2 + @parameters t, x, p[1:n_comp], q[1:n_comp] + @variables u(..)[1:n_comp] + Dt = Differential(t) + Dx = Differential(x) + Dxx = Differential(x)^2 + params = Symbolics.scalarize(reduce(vcat,[p .=> [1.5, 2.0], q .=> [1.2, 1.8]])) + # 1D PDE and boundary conditions + + eqs = [Dt(u(t, x)[i]) ~ p[i] * Dxx(u(t, x)[i]) for i in 1:n_comp] + + bcs = [[u(0, x)[i] ~ q[i] * cos(x), + u(t, 0)[i] ~ sin(t), + u(t, 1)[i] ~ exp(-t) * cos(1), + Dx(u(t,0)[i]) ~ 0.0] for i in 1:n_comp] + bcs_collected = reduce(vcat, bcs) + + # Space and time domains + domains = [t ∈ Interval(0.0, 1.0), + x ∈ Interval(0.0, 1.0)] + + # PDE system + + @named pdesys = PDESystem(eqs, bcs_collected, domains, [t, x], [u(t, x)[i] for i in 1:n_comp], Symbolics.scalarize(params)) + + + # Method of lines discretization + dx = 0.1 + order = 2 + discretization = MOLFiniteDifference([x => dx], t; approx_order = order) + + # Convert the PDE problem into an ODE problem + prob = discretize(pdesys,discretization) #error occurs here + + # Solve ODE problem + sol = solve(prob, Tsit5(), saveat=0.2) + + # Test that the system is correctly constructed + varname1 = Symbol("u_Any[1]") + varname2 = Symbol("u_Any[2]") + + + vars = @variables $varname1(..), $varname2(..) + + @test sol[u(t, x)[1]] == sol[vars[1](t, x)] + @test sol[u(t, x)[2]] == sol[vars[2](t, x)] + + @test sol(0.1, 0.1, dv = u(t, x)[1]) == sol(0.1, 0.1, dv = vars[1](t, x)) + @test sol(0.1, 0.1, dv = u(t, x)[2]) == sol(0.1, 0.1, dv = vars[2](t, x)) +end + +@testset "Budkyo-Sellers, nonlinlap case" begin + using MethodOfLines, DomainSets, ModelingToolkit, OrdinaryDiffEq + + @parameters φ t + @variables Tₛ(..) + + Dt = Differential(t) + Dφ = Differential(φ) + + A = 210 # W/m^2 + B = 2 # W/m^2K + D = 0.6 # W/m^2K + + P₂(x) = 0.5*(3*(x)^2 - 1) + α_func(ϕ) = 0.354 + 0.25*P₂(sin(ϕ)) + Q_func(t, ϕ) =430*cos(ϕ) + + s_in_y = 31556952.0 + + φmin = -Float64(π/2)+0.1 + φmax = Float64(π/2)-0.1 + + # Water fraction + f = 0.7 + # rho is density of Water + ρ = 1025 # kg/m^3 + # specific heat of Water + c_w = 4186 #J/(kgK) + # Depth of water + H = 70 #m + + #Heat Cap of earth + C = f*ρ*c_w*H + + eq = Dt(Tₛ(t, φ)) ~ ((1 - α_func(φ))*Q_func(t, φ) - A - B*Tₛ(t, φ) + D*Dφ(cos(φ)*Dφ(Tₛ(t, φ)))/cos(φ))/C + + bcs = [Tₛ(0, φ) ~ 12 - 40*P₂(sin(φ)), + Dφ(Tₛ(t, φmin)) ~ 0.0, + Dφ(Tₛ(t, φmax)) ~ 0.0] + domains = [t in IntervalDomain(0, 19*s_in_y), φ in IntervalDomain(φmin, φmax)] + + @named sys = PDESystem(eq, bcs, domains, [t, φ], [Tₛ(t, φ)]) + + discretization = MOLFiniteDifference([φ => 80], t, order = 4) + + prob = discretize(sys, discretization) # ERROR HERE + + sol = solve(prob, FBDF(), saveat = s_in_y) +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 53012a164..6d6f29c5a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,6 +21,7 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") if GROUP == "All" || GROUP == "MOL_Interface2" @time @safetestset "MOLFiniteDifference Interface" begin include("pde_systems/MOLtest2.jl") + include("pde_systems/MOLtest2_JSC.jl") end end @@ -107,6 +108,7 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") if GROUP == "All" || GROUP == "MOL_Interface1" @time @safetestset "MOLFiniteDifference Interface" begin include("pde_systems/MOLtest1.jl") + include("pde_systems/MOLtest1_JSC.jl") end end From a70f0d49774107754afe387bcd5c401c87d09c0f Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Mon, 4 Sep 2023 19:47:49 +0100 Subject: [PATCH 14/48] rebase (#4766) * add precompilation stage * merge precompiliation discretizations * fix * fix #293 * Update Project.toml * CompatHelper: add new compat entry for PrecompileTools at version 1, (keep existing compat) * Try new scheme for mixed derivs * remove show * mark broken * init mixed deriv rules * Update Project.toml --------- Co-authored-by: CompatHelper Julia --- Project.toml | 5 +- src/MethodOfLines.jl | 10 +- .../generate_finite_difference_rules.jl | 6 +- src/discretization/generate_ic_defaults.jl | 3 - .../2nd_order_mixed_deriv.jl | 27 ++++++ .../centered_difference.jl | 12 ++- src/precompile.jl | 81 ++++++++++++++++ .../pde_system_transformation.jl | 9 +- test/pde_systems/MOL_Mixed_Deriv.jl | 96 ++++++++++++++++++- test/runtests.jl | 11 +-- 10 files changed, 234 insertions(+), 26 deletions(-) create mode 100644 src/discretization/schemes/2nd_order_mixed_deriv/2nd_order_mixed_deriv.jl create mode 100644 src/precompile.jl diff --git a/Project.toml b/Project.toml index 79f8003f2..fb516f453 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MethodOfLines" uuid = "94925ecb-adb7-4558-8ed8-f975c56a0bf4" authors = ["Alex Jones, "] -version = "0.9.5" +version = "0.9.7" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" @@ -14,6 +14,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" PDEBase = "a7812802-0625-4b9e-961c-d332478797e5" +PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" @@ -35,6 +36,8 @@ Interpolations = "0.14" Latexify = "0.15, 0.16" ModelingToolkit = "8" OrdinaryDiffEq = "6" +PDEBase = "0.1.4" +PrecompileTools = "1" SafeTestsets = "0.0.1" SciMLBase = "1.94" PDEBase = "0.1.4" diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl index 5add161af..9065ae35e 100644 --- a/src/MethodOfLines.jl +++ b/src/MethodOfLines.jl @@ -12,11 +12,12 @@ using IfElse using StaticArrays using Interpolations using Latexify -import DomainSets +using PrecompileTools +using DomainSets # See here for the main `symbolic_discretize` and `generate_system` functions using PDEBase -using PDEBase: unitindices, unitindex, remove, insert, sym_dot, VariableMap, depvar, x2i, d_orders, vcat! +using PDEBase: unitindices, unitindex, remove, insert, sym_dot, VariableMap, depvar, x2i, d_orders, vcat!, update_varmap!, get_ops # To Extend import PDEBase.interface_errors import PDEBase.check_boundarymap @@ -78,6 +79,7 @@ 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/2nd_order_mixed_deriv/2nd_order_mixed_deriv.jl") include("discretization/schemes/upwind_difference/upwind_difference.jl") include("discretization/schemes/half_offset_centred_difference.jl") include("discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl") @@ -94,6 +96,10 @@ include("discretization/generate_ic_defaults.jl") include("scalar_discretization.jl") include("MOL_discretization.jl") +## PrecompileTools +include("precompile.jl") + +# Export export MOLFiniteDifference, discretize, symbolic_discretize, ODEFunctionExpr, generate_code, grid_align, edge_align, center_align, get_discrete, chebyspace export UpwindScheme, WENOScheme, FunctionalScheme diff --git a/src/discretization/generate_finite_difference_rules.jl b/src/discretization/generate_finite_difference_rules.jl index e443b83b3..8d2ad1b5d 100644 --- a/src/discretization/generate_finite_difference_rules.jl +++ b/src/discretization/generate_finite_difference_rules.jl @@ -50,7 +50,8 @@ function generate_finite_difference_rules(II::CartesianIndex, s::DiscreteSpace, if length(II) != 0 # Standard cartesian centered difference scheme central_deriv_rules_cartesian = generate_cartesian_rules(II, s, depvars, derivweights, bmap, indexmap, terms) - + # Mixed derivative rules + mixed_deriv_rules_cartesian = generate_mixed_rules(II, s, depvars, derivweights, bmap, indexmap, terms) # Advection rules if derivweights.advection_scheme isa UpwindScheme advection_rules = generate_winding_rules(II, s, depvars, derivweights, bmap, indexmap, terms) @@ -75,9 +76,10 @@ function generate_finite_difference_rules(II::CartesianIndex, s::DiscreteSpace, advection_rules = [] nonlinlap_rules = [] spherical_diffusion_rules = [] + mixed_deriv_rules_cartesian = [] integration_rules = [] end integration_rules = vcat(integration_rules, vec(generate_whole_domain_integration_rules(II, s, depvars, indexmap, terms))) - return vcat(vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules) + return vcat(vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(mixed_deriv_rules_cartesian), vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules) end diff --git a/src/discretization/generate_ic_defaults.jl b/src/discretization/generate_ic_defaults.jl index ecc24792f..dfc66fdc3 100644 --- a/src/discretization/generate_ic_defaults.jl +++ b/src/discretization/generate_ic_defaults.jl @@ -1,8 +1,5 @@ function PDEBase.generate_ic_defaults(tconds, s::DiscreteSpace, ::MOLFiniteDifference{G,DS}) where {G,DS<:ScalarizedDiscretization} t = s.time - if length(tconds) != length(s.ū) - @warn "The number of initial conditions must match the number of dependent variables." - end if s.time !== nothing u0 = mapreduce(vcat, tconds) do ic if isupper(ic) diff --git a/src/discretization/schemes/2nd_order_mixed_deriv/2nd_order_mixed_deriv.jl b/src/discretization/schemes/2nd_order_mixed_deriv/2nd_order_mixed_deriv.jl new file mode 100644 index 000000000..24303f569 --- /dev/null +++ b/src/discretization/schemes/2nd_order_mixed_deriv/2nd_order_mixed_deriv.jl @@ -0,0 +1,27 @@ +""" +Performs a mixed centered difference in `x` centered at index `II` of `u` +ufunc is a function that returns the correct discretization indexed at Itap, it is designed this way to allow for central differences of arbitrary expressions which may be needed in some schemes +""" +function mixed_central_difference((Dx, Dy), II, s, (xbs, ybs), (jx, ky), u, ufunc) + j, x = jx + k, y = ky + xweights, xItap = central_difference_weights_and_stencil(Dx, II, s, xbs, jx, u) + yweights, yItap = central_difference_weights_and_stencil(Dy, II, s, ybs, ky, u) + # TODO: Fix interface bcs + + out = sum(zip(xweights, xItap)) do (wx, xI) + sum(zip(yweights, yItap)) do (wy, yI) + xoffset = xI - II + yoffset = yI - II + I = II + xoffset + yoffset + wx * wy * ufunc(u, I, x) + end + end + + return out +end + +@inline function generate_mixed_rules(II::CartesianIndex, s::DiscreteSpace, depvars, derivweights::DifferentialDiscretizer, bcmap, indexmap, terms) + central_ufunc(u, I, x) = s.discvars[u][I] + return reduce(safe_vcat, [reduce(safe_vcat, [[(Differential(x)*Differential(y))(u) => mixed_central_difference((derivweights.map[Differential(x)], derivweights.map[Differential(y)]), Idx(II, s, u, indexmap), s, (filter_interfaces(bcmap[operation(u)][x]), filter_interfaces(bcmap[operation(u)][y])), ((x2i(s, u, x), x), (x2i(s, u, y), y)), u, central_ufunc) for y in remove(ivs(u, s), unwrap(x))] for x in ivs(u, s)], init = []) for u in depvars], init = []) +end \ No newline at end of file diff --git a/src/discretization/schemes/centered_difference/centered_difference.jl b/src/discretization/schemes/centered_difference/centered_difference.jl index b47d6dc62..64dc87903 100644 --- a/src/discretization/schemes/centered_difference/centered_difference.jl +++ b/src/discretization/schemes/centered_difference/centered_difference.jl @@ -2,7 +2,7 @@ Performs a centered difference in `x` centered at index `II` of `u` ufunc is a function that returns the correct discretization indexed at Itap, it is designed this way to allow for central differences of arbitrary expressions which may be needed in some schemes """ -function central_difference(D::DerivativeOperator{T,N,Wind,DX}, II, s, bs, jx, u, ufunc) where {T,N,Wind,DX<:Number} +function central_difference_weights_and_stencil(D::DerivativeOperator{T,N,Wind,DX}, II, s, bs, jx, u) where {T,N,Wind,DX<:Number} j, x = jx ndims(u, s) == 0 && return 0 # unit index in direction of the derivative @@ -23,10 +23,10 @@ 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 sym_dot(weights, ufunc(u, Itap, x)) + return weights, Itap end -function central_difference(D::DerivativeOperator{T,N,Wind,DX}, II, s, bs, jx, u, ufunc) where {T,N,Wind,DX<:AbstractVector} +function central_difference_weights_and_stencil(D::DerivativeOperator{T,N,Wind,DX}, II, s, bs, jx, u) where {T,N,Wind,DX<:AbstractVector} j, x = jx @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." ndims(u, s) == 0 && return 0 @@ -48,6 +48,12 @@ 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 weights, Itap +end + +function central_difference(D, II, s, bs, jx, u, ufunc) + j, x = jx + weights, Itap = central_difference_weights_and_stencil(D, II, s, bs, jx, u) return sym_dot(weights, ufunc(u, Itap, x)) end diff --git a/src/precompile.jl b/src/precompile.jl new file mode 100644 index 000000000..fc1d333b6 --- /dev/null +++ b/src/precompile.jl @@ -0,0 +1,81 @@ +@setup_workload begin + @compile_workload begin + begin + @parameters x y t + @variables u(..) v(..) + Dt = Differential(t) + Dx = Differential(x) + Dy = Differential(y) + Dxx = Differential(x)^2 + Dyy = Differential(y)^2 + + ∇²(u) = Dxx(u) + Dyy(u) + + brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.0 + + x_min = y_min = t_min = 0.0 + x_max = y_max = 1.0 + t_max = 5.0 + + α = 10.0 + + u0(x, y, t) = 22(y * (1 - y))^(3 / 2) + v0(x, y, t) = 27(x * (1 - x))^(3 / 2) + + eq = [Dt(u(x, y, t)) ~ 1.0 + v(x, y, t) * u(x, y, t)^2 - 4.4 * u(x, y, t) + α * ∇²(u(x, y, t)) + brusselator_f(x, y, t) + Dx(u(x, y, t)) + Dy(v(x, y, t)) + Dt(v(x, y, t)) ~ 3.4 * u(x, y, t) - v(x, y, t) * u(x, y, t)^2 + α * ∇²(v(x, y, t)) + Dx(v(x, y, t)*Dx(u(x, y, t))) + Dy(v(x, y, t)*Dy(u(x, y, t)))] + + domains = [x ∈ Interval(x_min, x_max), + y ∈ Interval(y_min, y_max), + t ∈ Interval(t_min, t_max)] + + # Periodic BCs + bcs = [u(x, y, 0) ~ u0(x, y, 0), + u(0, y, t) ~ u(1, y, t), + u(x, 0, t) ~ u(x, 1, t), v(x, y, 0) ~ v0(x, y, 0), + v(0, y, t) ~ v(1, y, t), + v(x, 0, t) ~ v(x, 1, t)] + + @named pdesys = PDESystem(eq, bcs, domains, [x, y, t], [u(x, y, t), v(x, y, t)]) + + N = 6 + + order = 2 + + discretization = MOLFiniteDifference([x => N, y => N], t, approx_order=order, grid_align=center_align) + + # Convert the PDE problem into an ODE problem + prob = discretize(pdesys, discretization) + + end + begin + @parameters x t + @variables u(..) + Dx = Differential(x) + Dt = Differential(t) + x_min = 0.0 + x_max = 1.0 + t_min = 0.0 + t_max = 6.0 + + analytic_u2(t, x) = x / (t + 1) + + eq = Dt(u(t, x)) ~ -u(t, x) * Dx(u(t, x)) + + bcs = [u(0, x) ~ x, + u(t, x_min) ~ analytic_u2(t, x_min), + u(t, x_max) ~ analytic_u2(t, x_max)] + + domains = [t ∈ Interval(t_min, t_max), + x ∈ Interval(x_min, x_max)] + + dx = 6 + + @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)]) + + disc = MOLFiniteDifference([x => dx], t, advection_scheme=WENOScheme(), grid_align=edge_align) + + prob = discretize(pdesys, disc) + end + end +end \ No newline at end of file diff --git a/src/system_parsing/pde_system_transformation.jl b/src/system_parsing/pde_system_transformation.jl index 1328ef22c..33f9ba7c8 100644 --- a/src/system_parsing/pde_system_transformation.jl +++ b/src/system_parsing/pde_system_transformation.jl @@ -5,7 +5,6 @@ Modified copilot explanation: """ function PDEBase.transform_pde_system!(v::PDEBase.VariableMap, boundarymap, sys::PDESystem, disc::MOLFiniteDifference) - eqs = copy(sys.eqs) bcs = copy(sys.bcs) done = false @@ -46,13 +45,15 @@ end """ Returns the term if it is incompatible, and whether to expand the term. """ -function filter_equivalent_differentials(term, differential, v) +function filter_differentials(term, differential, v, depth=0) S = Symbolics SU = SymbolicUtils if S.istree(term) op = SU.operation(term) if op isa Differential && isequal(op.x, differential.x) - return filter_equivalent_differentials(SU.arguments(term)[1], differential, v) + return filter_differentials(SU.arguments(term)[1], differential, v, depth + 1) + elseif op isa Differential && !isequal(op.x, differential.x) && depth <= 1 + return check_deriv_arg(arguments(term)[1], v) else return check_deriv_arg(term, v) end @@ -135,7 +136,7 @@ function descend_to_incompatible(term, v) if nonlinlapterm !== nothing badterm, shouldexpand = check_deriv_arg(nonlinlapterm, v) else - badterm, shouldexpand = filter_equivalent_differentials(term, op, v) + badterm, shouldexpand = filter_differentials(arguments(term)[1], op, v, 1) end if badterm !== nothing diff --git a/test/pde_systems/MOL_Mixed_Deriv.jl b/test/pde_systems/MOL_Mixed_Deriv.jl index 82c4b5261..cd67ae6e0 100644 --- a/test/pde_systems/MOL_Mixed_Deriv.jl +++ b/test/pde_systems/MOL_Mixed_Deriv.jl @@ -47,13 +47,13 @@ end @parameters t x y @variables u(..) Dt = Differential(t) - Dxy = Differential(x)*Differential(y) + Dxxy = Differential(x)^2 *Differential(y) - eq = [Dt(u(t, x, y)) ~ Dxy(u(t, x, y))] + eq = [Dt(u(t, x, y)) ~ Dxxy(u(t, x, y))] bcs = [u(0, x, y) ~ sinpi(x + y), - u(t, 0, y) ~ u(t, 1, y), - u(t, x, 0) ~ u(t, x, 1)] + u(t, 0, y) ~ sinpi(y), + u(t, x, 0) ~ sinpi(x)] domain = [t ∈ Interval(0.0, 1.0), x ∈ Interval(0.0, 1.0), @@ -67,5 +67,91 @@ end prob = discretize(pdesys, discretization, advection_scheme = WENOScheme()) sol = solve(prob, FBDF(), saveat = 0.1); - @test sol.retcode == SciMLBase.ReturnCode.Success + @test_broken sol.retcode == SciMLBase.ReturnCode.Success +end + +@test_broken begin#@testset "Mixed steady state problem" begin + @parameters x y + @variables u(..) + Dx = Differential(x) + Dy = Differential(y) + Dxx = Differential(x)^2 + Dyy = Differential(y)^2 + Dxy = Differential(x)*Differential(y) + + eq = [Dxx(u(x, y)) + Dyy(u(x, y)) + Dxy(u(x, y)) ~ 0] + + bcs = [u(0, y) ~ 0, + #Dx(u(0, y)) ~ y, + u(1, y) ~ y, + #Dx(u(1, y)) ~ y, + u(x, 0) ~ 0, + #Dy(u(x, 0)) ~ x, + u(x, 1) ~ x, + #Dy(u(x, 1)) ~ x + ] + + domain = [x ∈ Interval(0.0, 1.0), + y ∈ Interval(0.0, 1.0)] + + analytic_u(x, y) = x*y + + @named pdesys = PDESystem(eq, bcs, domain, [x, y], [u(x,y)]) + + dx = 0.1 + dy = 0.08 + + disc = MOLFiniteDifference([x => 20, y => 20], order = 4) + + prob = discretize(pdesys, disc) + + sol = solve(prob, NewtonRaphson()); + + @test_broken sol.retcode == SciMLBase.ReturnCode.Success + + solu = sol[u(x, y)] + solx = sol[x] + soly = sol[y] + + asol = [analytic_u(x, y) for x in solx[1:end-1], y in soly[1:end-1]] + + @test_broken solu[1:end-1, 1:end-1] ≈ asol atol = 1e-3 end + +@testset "Wave Equation u_tt ~ u_xx" begin + @parameters t x + @variables u(..) + Dt = Differential(t) + Dtt = Differential(t)^2 + Dxx = Differential(x)^2 + + eq = [Dtt(u(t, x)) ~ Dxx(u(t, x))] + + bcs = [u(0, x) ~ sinpi(x), + Dt(u(0, x)) ~ 0, + u(t, 0) ~ 0, + u(t, 1) ~ 0] + + domain = [t ∈ Interval(0.0, 1.0), + x ∈ Interval(0.0, 1.0)] + + analytic_u(t, x) = sinpi(x)*cospi(t) + + @named pdesys = PDESystem(eq, bcs, domain, [t, x], [u(t,x)]) + + dx = 0.01 + + disc = MOLFiniteDifference([x => dx], t) + + prob = discretize(pdesys, disc) + + sol = solve(prob, FBDF(), saveat = 0.1) + + xdisc = sol[x] + tdisc = sol[t] + usol = sol[u(t,x)] + + asol = [analytic_u(t, x) for t in tdisc, x in xdisc] + + @test usol ≈ asol atol = 1e-2 +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 6d6f29c5a..e8b83663b 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 == "Mixed_Derivatives" + @time @safetestset "MOLFiniteDifference Interface: Mixed Derivatives" begin + include("pde_systems/MOL_Mixed_Deriv.jl") + end + end if GROUP == "All" || GROUP == "Brusselator" @time @safetestset "MOLFiniteDifference Interface: 2D Brusselator Equation" begin include("pde_systems/brusselator_eq.jl") @@ -32,12 +37,6 @@ 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 == "Integrals" @time @safetestset "MOLFiniteDifference Interface: Integrals" begin include("pde_systems/MOL_1D_Integration.jl") From f44f965a6d900b2307e4ee155c9afcd00373ca54 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Wed, 6 Sep 2023 14:39:28 +0100 Subject: [PATCH 15/48] bump --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index fb516f453..bc9063210 100644 --- a/Project.toml +++ b/Project.toml @@ -40,7 +40,6 @@ PDEBase = "0.1.4" PrecompileTools = "1" SafeTestsets = "0.0.1" SciMLBase = "1.94" -PDEBase = "0.1.4" StaticArrays = "1" SymbolicUtils = "1" Symbolics = "5" From 941c298675cb7306d6f9360c9874d06b25092d3e Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Wed, 6 Sep 2023 15:21:10 +0100 Subject: [PATCH 16/48] Document --- docs/pages.jl | 1 + docs/src/performance.md | 10 ++++++++++ 2 files changed, 11 insertions(+) create mode 100644 docs/src/performance.md diff --git a/docs/pages.jl b/docs/pages.jl index fd89bdeab..8cb9687d6 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -9,6 +9,7 @@ pages = [ "tutorials/icbc_sampled.md", "tutorials/PIDE.md", ], + "Performance Tips" => "performance.md" "MOLFiniteDifference" => "MOLFiniteDifference.md", "Solution Interface - PDESolutions" => "solutions.md", "Grid and Solution Retrieval - Deprecated" => "get_grid.md", diff --git a/docs/src/performance.md b/docs/src/performance.md new file mode 100644 index 000000000..065148953 --- /dev/null +++ b/docs/src/performance.md @@ -0,0 +1,10 @@ +# Performance Tips + +While the default ModelingTolkit backend is highly capable, it doesn't scale well to very large systems of ODEs, present iun high resolution and/or higher dimensional PDE semidiscretizations that MethodOfLines.jl is capable of generating. + +To overcome this limitation, and scale MOL to realistic physics simulations, JuliaHub Inc. has released [JuliaSimCompiler.jl (installation instructions here)](https://help.juliahub.com/juliasimcompiler/dev/). + +Whenever this package and MethodOfLines are loaded together, MOL will also load `MethodOfLinesJuliaSimCompilerExt.jl`, an extension that takes advantage of this backend to speed up your compilation and solves automatically, with the same interface you are used to. + +!!! note + JuliaSimCompiler is part of JuliaSim and thus requires a valid JuliaSim license to use. JuliaSim is a proprietary software developed by JuliaHub Inc. Using the packages through the registry requires a valid JuliaSim license. It is free to use for non-commercial academic teaching and research purposes. For commercial users, license fees apply. Please refer to the [End User License Agreement](https://juliahub.com/company/eula/?_gl=1*120lqg6*_ga*MTAxODQ4OTE3Mi4xNjk0MDA4MDM5*_ga_8FC7JQQLXX*MTY5NDAwODAzOC4xLjEuMTY5NDAwODgxMC4wLjAuMA..) for details. Please contact [sales@juliahub.com](sales@juliahub.com) for purchasing information. \ No newline at end of file From 4245b4a6116aeac5fd7ccfba46dd69121c5d8646 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 11 Sep 2023 08:49:56 +0000 Subject: [PATCH 17/48] Bump actions/checkout from 3 to 4 Bumps [actions/checkout](https://github.com/actions/checkout) from 3 to 4. - [Release notes](https://github.com/actions/checkout/releases) - [Changelog](https://github.com/actions/checkout/blob/main/CHANGELOG.md) - [Commits](https://github.com/actions/checkout/compare/v3...v4) --- updated-dependencies: - dependency-name: actions/checkout dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] --- .github/workflows/CI.yml | 2 +- .github/workflows/Documentation.yml | 2 +- .github/workflows/Downstream.yml | 4 ++-- .github/workflows/Invalidations.yml | 4 ++-- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 9dd339e61..33cbab1ff 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -34,7 +34,7 @@ jobs: - '1' - '1.6' steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.version }} diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index d2b8c0fab..adc1628ef 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -11,7 +11,7 @@ jobs: build: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@latest with: version: '1' diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml index 9976af9c8..19cde18c4 100644 --- a/.github/workflows/Downstream.yml +++ b/.github/workflows/Downstream.yml @@ -19,14 +19,14 @@ jobs: package: - {user: SciML, repo: PDESystemLibrary.jl, group: MOL} steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.julia-version }} arch: x64 - uses: julia-actions/julia-buildpkg@latest - name: Clone Downstream - uses: actions/checkout@v3 + uses: actions/checkout@v4 with: repository: ${{ matrix.package.user }}/${{ matrix.package.repo }} path: downstream diff --git a/.github/workflows/Invalidations.yml b/.github/workflows/Invalidations.yml index 4d0004e83..28b9ce2fa 100644 --- a/.github/workflows/Invalidations.yml +++ b/.github/workflows/Invalidations.yml @@ -19,12 +19,12 @@ jobs: - uses: julia-actions/setup-julia@v1 with: version: '1' - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-invalidations@v1 id: invs_pr - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 with: ref: ${{ github.event.repository.default_branch }} - uses: julia-actions/julia-buildpkg@v1 From 770935ad36fa648c690d8bcae7dc26f469e7feab Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Mon, 11 Sep 2023 23:01:49 +0100 Subject: [PATCH 18/48] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 0ff9e0821..6b3d666f4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MethodOfLines" uuid = "94925ecb-adb7-4558-8ed8-f975c56a0bf4" authors = ["Alex Jones, "] -version = "0.9.7" +version = "0.9.6" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" From d62dcdf53a8886836f12cfde1031c9c92eabee91 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Fri, 15 Sep 2023 08:18:25 -0400 Subject: [PATCH 19/48] ci: add registry support for running JuliaSimCompiler integration tests Add a new workflow, `JuliaSimCompiler.yml`, which tests for integration with JuliaSimCompiler. Since PRs opened from forks don't have access to secrets, turn on `fail-fast` for these tests as well, so that it doesn't mess with the CI statuses for such PRs. --- .github/workflows/JuliaSimCompiler.yml | 52 +++++++++++++++++++ .../MethodOfLinesJuliaSimCompilerExt.jl | 4 +- src/interface/MOLFiniteDifference.jl | 5 +- test/pde_systems/MOLtest2_JSC.jl | 12 ++--- test/pde_systems/brusselator_eq.jl | 12 ++--- test/runtests.jl | 13 ++++- 6 files changed, 81 insertions(+), 17 deletions(-) create mode 100644 .github/workflows/JuliaSimCompiler.yml diff --git a/.github/workflows/JuliaSimCompiler.yml b/.github/workflows/JuliaSimCompiler.yml new file mode 100644 index 000000000..1cd21ad88 --- /dev/null +++ b/.github/workflows/JuliaSimCompiler.yml @@ -0,0 +1,52 @@ +name: JuliaSim Integration +on: + pull_request: + branches: + - master + push: + branches: + - master + tags: + - 'v*' +env: + JULIA_PKG_SERVER: https://juliahub.com/ +jobs: + test: + runs-on: ubuntu-latest + strategy: + fail-fast: true + matrix: + group: + - MOL_Interface1_JSC + - MOL_Interface2_JSC + version: + - '1' + - '1.6' + steps: + - uses: actions/checkout@v3 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + - uses: PumasAI/add-private-registry@main + with: + juliahub_token_encoded: ${{ secrets.JULIAHUB_TOKEN_ENCODED }} + private_registry_name: ${{ secrets.PRIVATE_REGISTRY_NAME }} + private_registry_uuid: ${{ secrets.PRIVATE_REGISTRY_UUID }} + - uses: actions/cache@v3 + env: + cache-name: cache-artifacts + with: + path: ~/.julia/artifacts + key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} + restore-keys: | + ${{ runner.os }}-test-${{ env.cache-name }}- + ${{ runner.os }}-test- + ${{ runner.os }}- + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 + env: + GROUP: ${{ matrix.group }} + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v3 + with: + file: lcov.info diff --git a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl index 0d58c23c0..03d2e008a 100644 --- a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl +++ b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl @@ -39,7 +39,9 @@ function PDEBase.generate_system(disc_state::PDEBase.EquationState, s, u0, tspan # * 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, checks = checks) - sys = IRSystem(sys) + if disc.useIR + sys = IRSystem(sys) + end return sys, tspan end catch e diff --git a/src/interface/MOLFiniteDifference.jl b/src/interface/MOLFiniteDifference.jl index b55c0b796..a4a7d7173 100644 --- a/src/interface/MOLFiniteDifference.jl +++ b/src/interface/MOLFiniteDifference.jl @@ -32,11 +32,12 @@ struct MOLFiniteDifference{G,D} <: AbstractEquationSystemDiscretization should_transform::Bool use_ODAE::Bool disc_strategy::D + useIR::Bool kwargs end # Constructors. If no order is specified, both upwind and centered differences will be 2nd order -function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, advection_scheme = UpwindScheme(), grid_align=CenterAlignedGrid(), discretization_strategy = ScalarizedDiscretization(), upwind_order = nothing, should_transform = true, use_ODAE = false, kwargs...) +function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, advection_scheme = UpwindScheme(), grid_align=CenterAlignedGrid(), discretization_strategy = ScalarizedDiscretization(), upwind_order = nothing, should_transform = true, use_ODAE = false, useIR = true, kwargs...) if upwind_order !== nothing @warn "`upwind_order` no longer does anything, and will be removed in a future release. See the docs for the current interface." end @@ -49,7 +50,7 @@ function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, advection_sche dxs = dxs isa Dict ? dxs : Dict(dxs) - return MOLFiniteDifference{typeof(grid_align), typeof(discretization_strategy)}(dxs, time, approx_order, advection_scheme, grid_align, should_transform, use_ODAE, discretization_strategy, kwargs) + return MOLFiniteDifference{typeof(grid_align), typeof(discretization_strategy)}(dxs, time, approx_order, advection_scheme, grid_align, should_transform, use_ODAE, discretization_strategy, useIR, kwargs) end PDEBase.get_time(disc::MOLFiniteDifference) = disc.time diff --git a/test/pde_systems/MOLtest2_JSC.jl b/test/pde_systems/MOLtest2_JSC.jl index 23c7e3f2c..dadf55315 100644 --- a/test/pde_systems/MOLtest2_JSC.jl +++ b/test/pde_systems/MOLtest2_JSC.jl @@ -120,7 +120,7 @@ using Test end -@testset "Nonlinlap with flux interface boundary condition" begin +@test_broken begin #@testset "Nonlinlap with flux interface boundary condition" begin @parameters t x1 x2 @variables c1(..) @variables c2(..) @@ -171,7 +171,7 @@ end @test sol.retcode == SciMLBase.ReturnCode.Success end -@testset "Another boundaries appearing in equations case" begin +@test_broken begin#@testset "Another boundaries appearing in equations case" begin g = 9.81 @@ -260,7 +260,7 @@ end end -@testset "Dt in BCs" begin +@test_broken begin #@testset "Dt in BCs" begin # Parameters, variables, and derivatives @parameters t x @variables u(..) @@ -296,7 +296,7 @@ end solu = sol[u(t, x)] # Temperature should increase with time end -@testset "ODE connected to PDE at boundary" begin +@test_broken begin #@testset "ODE connected to PDE at boundary" begin @variables u(..) v(..) w(..) @parameters t, r Dt = Differential(t) @@ -344,7 +344,7 @@ end solw = sol[w(t)] end -@testset "ODE connected to PDE" begin +@test_broken begin #@testset "ODE connected to PDE" begin @parameters t x @variables u(..) v(..) @@ -374,7 +374,7 @@ end solv = sol[v(t)] end -@testset "New style array variable conversion and interception" begin +@test_broken begin #@testset "New style array variable conversion and interception" begin # Parameters, variables, and derivatives n_comp = 2 @parameters t, x, p[1:n_comp], q[1:n_comp] diff --git a/test/pde_systems/brusselator_eq.jl b/test/pde_systems/brusselator_eq.jl index 263be3bb7..b1bd28e41 100644 --- a/test/pde_systems/brusselator_eq.jl +++ b/test/pde_systems/brusselator_eq.jl @@ -1,10 +1,10 @@ -using DomainSets, JuliaSimCompiler +using DomainSets using ModelingToolkit, MethodOfLines, LinearAlgebra, OrdinaryDiffEq # using Plots # local sol -begin #@testset "Test 01: Brusselator equation 2D" begin +@testset "Test 01: Brusselator equation 2D" begin @parameters x y t @variables u(..) v(..) Difft = Differential(t) @@ -61,6 +61,10 @@ begin #@testset "Test 01: Brusselator equation 2D" begin println("Solve:") @time sol = solve(prob, TRBDF2(),saveat=0.01) + solu = sol[u(x, y, t)] + solv = sol[v(x, y, t)] + + t = sol[t] # Solve reference problem xyd_brusselator = range(0,stop=1,length=N) @@ -97,10 +101,6 @@ begin #@testset "Test 01: Brusselator equation 2D" begin msol = solve(prob,TRBDF2(),saveat=0.01) # 2.771 s (5452 allocations: 65.73 MiB) - solu = sol[u(x, y, t)] - solv = sol[v(x, y, t)] - - t = sol[t] @testset "." begin for k in div(length(t), 2):length(t) msolu = msol.u[k][:,:,1] diff --git a/test/runtests.jl b/test/runtests.jl index 464be2922..25ea30719 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -27,7 +27,6 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") if GROUP == "All" || GROUP == "MOL_Interface2" @time @safetestset "MOLFiniteDifference Interface" begin include("pde_systems/MOLtest2.jl") - include("pde_systems/MOLtest2_JSC.jl") end end @@ -108,7 +107,6 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") if GROUP == "All" || GROUP == "MOL_Interface1" @time @safetestset "MOLFiniteDifference Interface" begin include("pde_systems/MOLtest1.jl") - include("pde_systems/MOLtest1_JSC.jl") end end @@ -129,4 +127,15 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") include("pde_systems/MOL_1D_Linear_Convection.jl") end end + ############### JSC ##################### + if GROUP == "All" || GROUP == "MOL_Interface2_JSC" + @time @safetestset "MOLFiniteDifference Interface" begin + include("pde_systems/MOLtest2_JSC.jl") + end + end + if GROUP == "All" || GROUP == "MOL_Interface1_JSC" + @time @safetestset "MOLFiniteDifference Interface" begin + include("pde_systems/MOLtest1_JSC.jl") + end + end end From 54e236e66fd3788b11f1f15a2da502ce5324f2a1 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Mon, 18 Sep 2023 14:13:09 +0100 Subject: [PATCH 20/48] add option --- src/MOL_discretization.jl | 2 +- src/MethodOfLines.jl | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/MOL_discretization.jl b/src/MOL_discretization.jl index 72caa5d52..26bf1dd98 100644 --- a/src/MOL_discretization.jl +++ b/src/MOL_discretization.jl @@ -97,7 +97,7 @@ function SciMLBase.ODEFunction(pdesys::PDESystem, discretization::MethodOfLines. gridlocs = get_gridloc.(us, (s,)) f_analytic = generate_function_from_gridlocs(analytic, gridlocs, s) end - return ODEFunction(simpsys; analytic = f_analytic, discretization.kwargs..., kwargs...) + return ODEFunction(simpsys; analytic = f_analytic, eval_module = @__MODULE__, discretization.kwargs..., kwargs...) end catch e println("The system of equations is:") diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl index 9065ae35e..7a818230e 100644 --- a/src/MethodOfLines.jl +++ b/src/MethodOfLines.jl @@ -14,6 +14,9 @@ using Interpolations using Latexify using PrecompileTools using DomainSets +using RuntimeGeneratedFunctions +RuntimeGeneratedFunctions.init(@__MODULE__) + # See here for the main `symbolic_discretize` and `generate_system` functions using PDEBase From 278a5012af048bd30cb6f22f6a18769a4fd84181 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Mon, 18 Sep 2023 16:22:20 +0100 Subject: [PATCH 21/48] mark broken --- test/pde_systems/MOLtest1_JSC.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/pde_systems/MOLtest1_JSC.jl b/test/pde_systems/MOLtest1_JSC.jl index 6e31c3a71..0f4d13a3b 100644 --- a/test/pde_systems/MOLtest1_JSC.jl +++ b/test/pde_systems/MOLtest1_JSC.jl @@ -29,7 +29,7 @@ using Test sol = solve(prob, Tsit5()) end -@testset "Heat Equation 2D 1 variable" begin +@test_broken begin #@testset "Heat Equation 2D 1 variable" begin @parameters t x y @variables u(..) Dxx = Differential(x)^2 @@ -316,7 +316,7 @@ end end -@testset "2D variable connected to 1D variable at boundary #33" begin +@test_broken begin #@testset "2D variable connected to 1D variable at boundary #33" begin @parameters t x r @variables u(..) v(..) Dt = Differential(t) From b64af4001574279d4ce06364b5d63861d6f39c14 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Mon, 18 Sep 2023 16:23:20 +0100 Subject: [PATCH 22/48] add rgf --- Project.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index bc9063210..bff87bfc4 100644 --- a/Project.toml +++ b/Project.toml @@ -15,6 +15,7 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" PDEBase = "a7812802-0625-4b9e-961c-d332478797e5" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" @@ -47,10 +48,10 @@ TermInterface = "0.2, 0.3" julia = "1.6" [extras] +JuliaSimCompiler = "8391cb6b-4921-5777-4e45-fd9aab8cb88d" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" -JuliaSimCompiler = "8391cb6b-4921-5777-4e45-fd9aab8cb88d" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] From d88bf45c7692c23a41536b6a9a49ae8286f813b1 Mon Sep 17 00:00:00 2001 From: Anant Thazhemadam Date: Wed, 20 Sep 2023 09:07:48 +0530 Subject: [PATCH 23/48] test: add JuliaSimCompiler only to run its integration tests Remove `JuliaSimCompiler` as a test dependency as that would require having access to the package at all times to run the tests, making it synonymous with a direct dependency (and less of an extension) in the context of running tests. Instead, explicitly `Pkg.add` `JuliaSimCompiler` instead in the relevant test matrices. --- Project.toml | 4 ++-- test/runtests.jl | 3 +++ 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index bff87bfc4..2d35eb7a8 100644 --- a/Project.toml +++ b/Project.toml @@ -48,11 +48,11 @@ TermInterface = "0.2, 0.3" julia = "1.6" [extras] -JuliaSimCompiler = "8391cb6b-4921-5777-4e45-fd9aab8cb88d" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "NonlinearSolve", "SafeTestsets", "StableRNGs", "JuliaSimCompiler"] +test = ["Test", "NonlinearSolve", "SafeTestsets", "StableRNGs", "Pkg"] diff --git a/test/runtests.jl b/test/runtests.jl index 25ea30719..7c3e8be40 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,5 @@ using SafeTestsets +using Pkg const GROUP = get(ENV, "GROUP", "All") const is_APPVEYOR = Sys.iswindows() && haskey(ENV, "APPVEYOR") @@ -129,11 +130,13 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") end ############### JSC ##################### if GROUP == "All" || GROUP == "MOL_Interface2_JSC" + Pkg.add("JuliaSimCompiler") @time @safetestset "MOLFiniteDifference Interface" begin include("pde_systems/MOLtest2_JSC.jl") end end if GROUP == "All" || GROUP == "MOL_Interface1_JSC" + Pkg.add("JuliaSimCompiler") @time @safetestset "MOLFiniteDifference Interface" begin include("pde_systems/MOLtest1_JSC.jl") end From d7711464f4809750fa66bc1c3b0aba7dd2e7956c Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Wed, 20 Sep 2023 20:42:16 +0100 Subject: [PATCH 24/48] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 938a17da9..304b430d6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MethodOfLines" uuid = "94925ecb-adb7-4558-8ed8-f975c56a0bf4" authors = ["Alex Jones, "] -version = "0.9.6" +version = "0.10.0" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" From c32dcbb62ae3539a70e02c6e123c9af2acb8736e Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Fri, 22 Sep 2023 02:48:24 +0000 Subject: [PATCH 25/48] CompatHelper: bump compat for SciMLBase to 2, (keep existing compat) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 304b430d6..5e5b1f29e 100644 --- a/Project.toml +++ b/Project.toml @@ -40,7 +40,7 @@ OrdinaryDiffEq = "6" PDEBase = "0.1.4" PrecompileTools = "1" SafeTestsets = "0.0.1" -SciMLBase = "1.94" +SciMLBase = "1.94, 2" StaticArrays = "1" SymbolicUtils = "1" Symbolics = "5" From 7ca9bac8c42f2f7b5ad8fb27f6ef1db78493e18c Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 22 Sep 2023 00:10:45 -0400 Subject: [PATCH 26/48] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 5e5b1f29e..7a6dd3cbe 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MethodOfLines" uuid = "94925ecb-adb7-4558-8ed8-f975c56a0bf4" authors = ["Alex Jones, "] -version = "0.10.0" +version = "0.10.1" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" From ae411ea673f2ca39278a1f2b1383cb1bb4c4dac7 Mon Sep 17 00:00:00 2001 From: Arno Strouwen Date: Mon, 25 Sep 2023 03:59:02 +0200 Subject: [PATCH 27/48] Documenter 1.0 upgrade --- .github/workflows/CI.yml | 4 ++++ docs/Project.toml | 3 +-- docs/make.jl | 15 +++------------ docs/pages.jl | 2 +- docs/src/howitworks.md | 2 +- docs/src/index.md | 33 +++++++++++++-------------------- docs/src/performance.md | 2 +- 7 files changed, 24 insertions(+), 37 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 33cbab1ff..f87b2bd42 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -3,9 +3,13 @@ on: pull_request: branches: - master + paths-ignore: + - 'docs/**' push: branches: - master + paths-ignore: + - 'docs/**' jobs: test: runs-on: ubuntu-latest diff --git a/docs/Project.toml b/docs/Project.toml index ca8b8f6f8..1c1853490 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -10,9 +10,8 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" [compat] DifferentialEquations = "7.6" -Documenter = "0.27" +Documenter = "1" DomainSets = "0.5, 0.6" -MethodOfLines = "0.8, 0.9" ModelingToolkit = "8" NonlinearSolve = "1" OrdinaryDiffEq = "6" diff --git a/docs/make.jl b/docs/make.jl index be9a3fa46..63b7e3d9e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -10,19 +10,10 @@ include("pages.jl") makedocs(sitename = "MethodOfLines.jl", authors = "Chris Rackauckas, Alex Jones et al.", - clean = true, - doctest = false, - strict = [ - :doctest, - :linkcheck, - :parse_error, - :example_block, - # Other available options are - # :autodocs_block, :cross_references, :docs_block, :eval_block, :example_block, :footnote, :meta_block, :missing_docs, :setup_block, - ], + clean = true, doctest = false, linkcheck = true, modules = [MethodOfLines], - format = Documenter.HTML(analytics = "UA-90474609-3", - assets = ["assets/favicon.ico"], + warnonly = [:docs_block, :missing_docs, :cross_references], + format = Documenter.HTML(assets = ["assets/favicon.ico"], canonical = "https://docs.sciml.ai/MethodOfLines/stable/"), pages = pages) diff --git a/docs/pages.jl b/docs/pages.jl index 8cb9687d6..11caa32a2 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -9,7 +9,7 @@ pages = [ "tutorials/icbc_sampled.md", "tutorials/PIDE.md", ], - "Performance Tips" => "performance.md" + "Performance Tips" => "performance.md", "MOLFiniteDifference" => "MOLFiniteDifference.md", "Solution Interface - PDESolutions" => "solutions.md", "Grid and Solution Retrieval - Deprecated" => "get_grid.md", diff --git a/docs/src/howitworks.md b/docs/src/howitworks.md index 2557af4ad..db213afaf 100644 --- a/docs/src/howitworks.md +++ b/docs/src/howitworks.md @@ -6,7 +6,7 @@ See [here](https://github.com/SciML/MethodOfLines.jl/blob/master/src/MOL_discret Given your discretization and `PDESystem`, we take each independent variable defined on the space to be discretized and create a corresponding range. We then take each dependent variable and create an array of symbolic variables to represent it in its discretized form. This is stored in a [`DiscreteSpace`](https://github.com/SciML/MethodOfLines.jl/blob/master/src/discretization/discretize_vars.jl) object, a useful abstraction. -We recognize boundary conditions, i.e. whether they are on the upper or lower ends of the domain, or periodic [here](https://github.com/SciML/MethodOfLines.jl/blob/master/src/system_parsing/bcs/parse_boundaries.jl), and use this information to construct the interior of the domain for each equation [here](https://github.com/SciML/MethodOfLines.jl/blob/master/src/system_parsing/interiormap.jl). Each PDE is matched to each dependent variable in this step by which variable is of highest order in each PDE, with precedence given to time derivatives. This dictates which boundary conditions reduce the size of the interior for which PDE. This is done to ensure that there will be the same number of equations as discrete variable states, so that the system of equations is balanced. +We recognize boundary conditions, i.e. whether they are on the upper or lower ends of the domain, or periodic [here](https://github.com/SciML/PDEBase.jl/blob/master/src/parse_boundaries.jl), and use this information to construct the interior of the domain for each equation [here](https://github.com/SciML/MethodOfLines.jl/blob/master/src/system_parsing/interior_map.jl). Each PDE is matched to each dependent variable in this step by which variable is of highest order in each PDE, with precedence given to time derivatives. This dictates which boundary conditions reduce the size of the interior for which PDE. This is done to ensure that there will be the same number of equations as discrete variable states, so that the system of equations is balanced. Next, the boundary conditions are discretized, creating an equation for each point on the boundary in terms of the discretized variables, replacing any space derivatives in the direction of the boundary with their upwind finite difference expressions. [This](https://github.com/SciML/MethodOfLines.jl/blob/master/src/discretization/generate_bc_eqs.jl) is the place to look to see how this happens. diff --git a/docs/src/index.md b/docs/src/index.md index 58c7dba1c..80c1ffc40 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -109,26 +109,19 @@ Pkg.status(;mode = PKGMODE_MANIFEST) # hide ```@raw html ``` -```@raw html -You can also download the -manifest file and the -project file. +using Markdown +version = TOML.parse(read("../../Project.toml", String))["version"] +name = TOML.parse(read("../../Project.toml", String))["name"] +link_manifest = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version * + "/assets/Manifest.toml" +link_project = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version * + "/assets/Project.toml" +Markdown.parse("""You can also download the +[manifest]($link_manifest) +file and the +[project]($link_project) +file. +""") ``` diff --git a/docs/src/performance.md b/docs/src/performance.md index 065148953..77bd7e9b7 100644 --- a/docs/src/performance.md +++ b/docs/src/performance.md @@ -7,4 +7,4 @@ To overcome this limitation, and scale MOL to realistic physics simulations, Jul Whenever this package and MethodOfLines are loaded together, MOL will also load `MethodOfLinesJuliaSimCompilerExt.jl`, an extension that takes advantage of this backend to speed up your compilation and solves automatically, with the same interface you are used to. !!! note - JuliaSimCompiler is part of JuliaSim and thus requires a valid JuliaSim license to use. JuliaSim is a proprietary software developed by JuliaHub Inc. Using the packages through the registry requires a valid JuliaSim license. It is free to use for non-commercial academic teaching and research purposes. For commercial users, license fees apply. Please refer to the [End User License Agreement](https://juliahub.com/company/eula/?_gl=1*120lqg6*_ga*MTAxODQ4OTE3Mi4xNjk0MDA4MDM5*_ga_8FC7JQQLXX*MTY5NDAwODAzOC4xLjEuMTY5NDAwODgxMC4wLjAuMA..) for details. Please contact [sales@juliahub.com](sales@juliahub.com) for purchasing information. \ No newline at end of file + JuliaSimCompiler is part of JuliaSim and thus requires a valid JuliaSim license to use. JuliaSim is a proprietary software developed by JuliaHub Inc. Using the packages through the registry requires a valid JuliaSim license. It is free to use for non-commercial academic teaching and research purposes. For commercial users, license fees apply. Please refer to the [End User License Agreement](https://juliahub.com/company/eula/?_gl=1*120lqg6*_ga*MTAxODQ4OTE3Mi4xNjk0MDA4MDM5*_ga_8FC7JQQLXX*MTY5NDAwODAzOC4xLjEuMTY5NDAwODgxMC4wLjAuMA..) for details. Please contact [sales@juliahub.com](https://juliahub.com/products/pricing/) for purchasing information. \ No newline at end of file From e32aefd91b6c8a57abcc7434d9861d9fa9a719e1 Mon Sep 17 00:00:00 2001 From: Arno Strouwen Date: Mon, 25 Sep 2023 04:00:11 +0200 Subject: [PATCH 28/48] add formatting ci workflow --- .JuliaFormatter.toml | 2 ++ .github/workflows/FormatCheck.yml | 42 +++++++++++++++++++++++++++++++ 2 files changed, 44 insertions(+) create mode 100644 .JuliaFormatter.toml create mode 100644 .github/workflows/FormatCheck.yml diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 000000000..96727c8f1 --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1,2 @@ +style = "sciml" +format_markdown = true diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml new file mode 100644 index 000000000..dd551501c --- /dev/null +++ b/.github/workflows/FormatCheck.yml @@ -0,0 +1,42 @@ +name: format-check + +on: + push: + branches: + - 'master' + - 'release-' + tags: '*' + pull_request: + +jobs: + build: + runs-on: ${{ matrix.os }} + strategy: + matrix: + julia-version: [1] + julia-arch: [x86] + os: [ubuntu-latest] + steps: + - uses: julia-actions/setup-julia@latest + with: + version: ${{ matrix.julia-version }} + + - uses: actions/checkout@v4 + - name: Install JuliaFormatter and format + # This will use the latest version by default but you can set the version like so: + # + # julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="0.13.0"))' + run: | + julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))' + julia -e 'using JuliaFormatter; format(".", verbose=true)' + - name: Format check + run: | + julia -e ' + out = Cmd(`git diff --name-only`) |> read |> String + if out == "" + exit(0) + else + @error "Some files have not been formatted !!!" + write(stdout, out) + exit(1) + end' From 36ef0fd2cf79379c16cf11fb2c8b71f4243139a5 Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Mon, 25 Sep 2023 03:09:31 +0000 Subject: [PATCH 29/48] CompatHelper: add new compat entry for MethodOfLines at version 0.9 for package docs, (keep existing compat) --- docs/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/Project.toml b/docs/Project.toml index 1c1853490..e1642e490 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -12,6 +12,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" DifferentialEquations = "7.6" Documenter = "1" DomainSets = "0.5, 0.6" +MethodOfLines = "0.9" ModelingToolkit = "8" NonlinearSolve = "1" OrdinaryDiffEq = "6" From 5513b9f3d76d11a6424aafa3878c7e0f87b043ca Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 25 Sep 2023 08:18:42 +0000 Subject: [PATCH 30/48] Bump actions/checkout from 3 to 4 Bumps [actions/checkout](https://github.com/actions/checkout) from 3 to 4. - [Release notes](https://github.com/actions/checkout/releases) - [Changelog](https://github.com/actions/checkout/blob/main/CHANGELOG.md) - [Commits](https://github.com/actions/checkout/compare/v3...v4) --- updated-dependencies: - dependency-name: actions/checkout dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] --- .github/workflows/JuliaSimCompiler.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/JuliaSimCompiler.yml b/.github/workflows/JuliaSimCompiler.yml index 1cd21ad88..d042a1f1b 100644 --- a/.github/workflows/JuliaSimCompiler.yml +++ b/.github/workflows/JuliaSimCompiler.yml @@ -23,7 +23,7 @@ jobs: - '1' - '1.6' steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.version }} From 186809073a576797944dad92e30cdf93450af699 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Thu, 21 Sep 2023 17:18:38 +0100 Subject: [PATCH 31/48] move extension over --- Project.toml | 6 -- .../MethodOfLinesJuliaSimCompilerExt.jl | 90 ------------------- 2 files changed, 96 deletions(-) delete mode 100644 ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl diff --git a/Project.toml b/Project.toml index 7a6dd3cbe..affa1b02b 100644 --- a/Project.toml +++ b/Project.toml @@ -22,12 +22,6 @@ SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" TermInterface = "8ea1fca8-c5ef-4a55-8b96-4e9afe9c9a3c" -[weakdeps] -JuliaSimCompiler = "8391cb6b-4921-5777-4e45-fd9aab8cb88d" - -[extensions] -MethodOfLinesJuliaSimCompilerExt = ["JuliaSimCompiler"] - [compat] Combinatorics = "1" DiffEqBase = "6" diff --git a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl b/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl deleted file mode 100644 index 03d2e008a..000000000 --- a/ext/MethodOfLinesJuliaSimCompilerExt/MethodOfLinesJuliaSimCompilerExt.jl +++ /dev/null @@ -1,90 +0,0 @@ -module MethodOfLinesJuliaSimCompilerExt -using JuliaSimCompiler, MethodOfLines, PDEBase, SciMLBase, ModelingToolkit -using PDEBase: get_discvars, get_time, add_metadata! - -function PDEBase.add_metadata!(sys::JuliaSimCompiler.IRSystem, meta) - sys.info.parent.metadata.metadata[] = meta -end - -function PDEBase.generate_system(disc_state::PDEBase.EquationState, s, u0, tspan, metadata, - disc::MethodOfLines.MOLFiniteDifference) - println("JuliaSimCompiler: generate_system") - discvars = get_discvars(s) - t = get_time(disc) - name = metadata.pdesys.name - pdesys = metadata.pdesys - alleqs = vcat(disc_state.eqs, unique(disc_state.bceqs)) - alldepvarsdisc = vec(reduce(vcat, vec(unique(reduce(vcat, vec.(values(discvars))))))) - - defaults = Dict(pdesys.ps === nothing || pdesys.ps === SciMLBase.NullParameters() ? u0 : - vcat(u0, pdesys.ps)) - ps = pdesys.ps === nothing || pdesys.ps === SciMLBase.NullParameters() ? Num[] : - first.(pdesys.ps) - # Finalize - # 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: - # 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, 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, checks = checks) - if disc.useIR - sys = IRSystem(sys) - end - return sys, tspan - end - catch e - println("The system of equations is:") - println(alleqs) - println() - println("Discretization failed, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.") - println() - rethrow(e) - end -end - -function SciMLBase.discretize(pdesys::PDESystem, - discretization::MOLFiniteDifference; - analytic = nothing, kwargs...) - sys, tspan = SciMLBase.symbolic_discretize(pdesys, discretization) - try - simpsys = structural_simplify(sys) - if tspan === nothing - return prob = NonlinearProblem(simpsys, ones(length(simpsys.states)); - discretization.kwargs..., kwargs...) - else - # Use ODAE if nessesary - - prob = ODEProblem(simpsys, Pair[], tspan; discretization.kwargs..., - kwargs...) - if analytic === nothing - return prob - else - f = ODEFunction(pdesys, discretization, analytic = analytic, - discretization.kwargs..., kwargs...) - - return ODEProblem(f, prob.u0, prob.tspan, prob.p; - discretization.kwargs..., kwargs...) - end - end - catch e - PDEBase.error_analysis(sys, e) - end -end - -function PDEBase.error_analysis(sys::IRSystem, e) - rethrow(e) - -end -end \ No newline at end of file From 685bfd44a4c7dbd653ef51e0bfd33a935bdc2ee0 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Thu, 21 Sep 2023 17:39:52 +0100 Subject: [PATCH 32/48] remove tests --- .github/workflows/JuliaSimCompiler.yml | 52 --- test/pde_systems/MOLtest1_JSC.jl | 353 ------------------ test/pde_systems/MOLtest2_JSC.jl | 478 ------------------------- test/runtests.jl | 13 - 4 files changed, 896 deletions(-) delete mode 100644 .github/workflows/JuliaSimCompiler.yml delete mode 100644 test/pde_systems/MOLtest1_JSC.jl delete mode 100644 test/pde_systems/MOLtest2_JSC.jl diff --git a/.github/workflows/JuliaSimCompiler.yml b/.github/workflows/JuliaSimCompiler.yml deleted file mode 100644 index d042a1f1b..000000000 --- a/.github/workflows/JuliaSimCompiler.yml +++ /dev/null @@ -1,52 +0,0 @@ -name: JuliaSim Integration -on: - pull_request: - branches: - - master - push: - branches: - - master - tags: - - 'v*' -env: - JULIA_PKG_SERVER: https://juliahub.com/ -jobs: - test: - runs-on: ubuntu-latest - strategy: - fail-fast: true - matrix: - group: - - MOL_Interface1_JSC - - MOL_Interface2_JSC - version: - - '1' - - '1.6' - steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 - with: - version: ${{ matrix.version }} - - uses: PumasAI/add-private-registry@main - with: - juliahub_token_encoded: ${{ secrets.JULIAHUB_TOKEN_ENCODED }} - private_registry_name: ${{ secrets.PRIVATE_REGISTRY_NAME }} - private_registry_uuid: ${{ secrets.PRIVATE_REGISTRY_UUID }} - - uses: actions/cache@v3 - env: - cache-name: cache-artifacts - with: - path: ~/.julia/artifacts - key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} - restore-keys: | - ${{ runner.os }}-test-${{ env.cache-name }}- - ${{ runner.os }}-test- - ${{ runner.os }}- - - uses: julia-actions/julia-buildpkg@v1 - - uses: julia-actions/julia-runtest@v1 - env: - GROUP: ${{ matrix.group }} - - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v3 - with: - file: lcov.info diff --git a/test/pde_systems/MOLtest1_JSC.jl b/test/pde_systems/MOLtest1_JSC.jl deleted file mode 100644 index 0f4d13a3b..000000000 --- a/test/pde_systems/MOLtest1_JSC.jl +++ /dev/null @@ -1,353 +0,0 @@ -using ModelingToolkit, MethodOfLines, LinearAlgebra, OrdinaryDiffEq -using ModelingToolkit: operation, istree, arguments -using DomainSets -using NonlinearSolve -using StableRNGs -using JuliaSimCompiler -using Test - -# # Define some variables -@testset "Heat Equation 1D 2 variables" begin - @parameters t x - @variables u(..) v(..) - Dt = Differential(t) - Dx = Differential(x) - Dxx = Differential(x)^2 - eqs = [Dt(u(t, x)) ~ Dxx(u(t, x)), - Dt(v(t, x)) ~ Dxx(v(t, x))] - bcs = [u(0, x) ~ -x * (x - 1) * sin(x), - v(0, x) ~ -x * (x - 1) * sin(x), - u(t, 0) ~ 0.0, u(t, 1) ~ 0.0, - v(t, 0) ~ 0.0, v(t, 1) ~ 0.0] - - domains = [t ∈ Interval(0.0, 1.0), - x ∈ Interval(0.0, 1.0)] - - @named pdesys = PDESystem(eqs, bcs, domains, [t, x], [u(t, x), v(t, x)]) - discretization = MOLFiniteDifference([x => 0.1], t; grid_align=edge_align) - prob = discretize(pdesys, discretization) # This gives an ODEProblem since it's time-dependent - sol = solve(prob, Tsit5()) -end - -@test_broken begin #@testset "Heat Equation 2D 1 variable" begin - @parameters t x y - @variables u(..) - Dxx = Differential(x)^2 - Dyy = Differential(y)^2 - Dt = Differential(t) - t_min = 0.0 - t_max = 2.0 - x_min = 0.0 - x_max = 2.0 - y_min = 0.0 - y_max = 2.0 - - # 3D PDE - eq = Dt(u(t, x, y)) ~ Dxx(u(t, x, y)) + Dyy(u(t, x, y)) - - analytic_sol_func(t, x, y) = exp(x + y) * cos(x + y + 4t) - # Initial and boundary conditions - bcs = [u(t_min, x, y) ~ analytic_sol_func(t_min, x, y), - u(t, x_min, y) ~ analytic_sol_func(t, x_min, y), - u(t, x_max, y) ~ analytic_sol_func(t, x_max, y), - u(t, x, y_min) ~ analytic_sol_func(t, x, y_min), - u(t, x, y_max) ~ analytic_sol_func(t, x, y_max)] - - # Space and time domains - domains = [t ∈ Interval(t_min, t_max), - x ∈ Interval(x_min, x_max), - y ∈ Interval(y_min, y_max)] - @named pdesys = PDESystem([eq], bcs, domains, [t, x, y], [u(t, x, y)]) - # Method of lines discretization - dx = 0.1 - dy = 0.2 - discretization = MOLFiniteDifference([x => dx, y => dy], t) - prob = ModelingToolkit.discretize(pdesys, discretization) - sol = solve(prob, Tsit5()) -end - -# Diffusion in a sphere -@testset "Spherical Diffusion" begin - - @parameters t r - @variables u(..) - Dt = Differential(t) - Dr = Differential(r) - Drr = Dr^2 - eq = Dt(u(t, r)) ~ 1 / r^2 * Dr(r^2 * Dr(u(t, r))) - bcs = [u(0, r) ~ -r * (r - 1) * sin(r), - Dr(u(t, 0)) ~ 0.0, u(t, 1) ~ sin(1)] - - domains = [t ∈ Interval(0.0, 1.0), - r ∈ Interval(0.0, 1.0)] - - @named pdesys = PDESystem(eq, bcs, domains, [t, r], [u(t, r)]) - discretization = MOLFiniteDifference([r => 0.1], t) - prob = discretize(pdesys, discretization) # This gives an ODEProblem since it's time-dependent - sol = solve(prob, Tsit5()) -end - -@testset "RHS = 0" begin - @parameters x, t - @variables u(..) - Dt = Differential(t) - Dx = Differential(x) - Dx2 = Differential(x)^2 - Dx3 = Differential(x)^3 - Dx4 = Differential(x)^4 - - α = 1 - β = 4 - γ = 1 - eq = Dt(u(x, t)) + u(x, t) * Dx(u(x, t)) + α * Dx2(u(x, t)) + β * Dx3(u(x, t)) + γ * Dx4(u(x, t)) ~ 0 - - du(x, t; z=-x / 2 + t) = 15 / 2 * (tanh(z) + 1) * (3 * tanh(z) - 1) * sech(z)^2 - - bcs = [u(x, 0) ~ x^2, - Dx(u(-10, t)) ~ du(-10, t), - Dx(u(10, t)) ~ du(10, t)] - - # Space and time domains - domains = [x ∈ Interval(-10.0, 10.0), - t ∈ Interval(0.0, 0.5)] - # Discretization - dx = 0.4 - dt = 0.2 - - discretization = MOLFiniteDifference([x => dx], t; approx_order=4) - @named pdesys = PDESystem(eq, bcs, domains, [x, t], [u(x, t)]) - prob = discretize(pdesys, discretization) -end - -@testset "Wave Equation" begin - # Parameters, variables, and derivatives - @parameters t x - @variables u(..) - Dt = Differential(t) # Required in ICs and equation - Dtt = Differential(t)^2 # required in equation - Dxxxx = Differential(x)^4 # required in equation - Dxx = Differential(x)^2 # required in BCs - # some parameters - EI = 291.6667 - m = 1.3850 - c = 1.0 - p = 0.01 - L = 2.0 - - # 1D PDE and boundaru conditions - eq = m * Dtt(u(t, x)) + c * Dt(u(t, x)) + EI * Dxxxx(u(t, x)) ~ 0 - - ic_bc = [u(0, x) ~ (p * x * (x^3 + L^3 - 2 * L * x^2) / (24 * EI)), #for all 0 < u < L - Dt(u(0, x)) ~ 0.0, # for all 0 < u < L - u(t, 0) ~ 0.0, # for all t > 0,, displacement zero at u=0 - u(t, 2) ~ 0.0, # for all t > 0,, displacement zero at u=L - Dxx(u(t, 0)) ~ 0.0, # for all t > 0,, curvature zero at u=0 - Dxx(u(t, 2)) ~ 0.0] # for all t > 0,, curvature zero at u=L - - # Space and time domains - domains = [t ∈ Interval(0.0, 2.0), - x ∈ Interval(0.0, L)] - - dt = 0.1 # dt related to saving the data.. not actual dt - # PDE sustem - @named pdesys = PDESystem(eq, ic_bc, domains, [t, x], [u(t, x)]) - - # Method of lines discretization - dx = 0.1 - order = 2 - discretization = MOLFiniteDifference([x => dx, t => dt], approx_order=order) - - # Convert the PDE problem into an ODE problem - prob = discretize(pdesys, discretization) - - # Solve the ODE problem - sol = NonlinearSolve.solve(prob, NewtonRaphson()) -end - -@testset "Rearranged Robin" begin - @parameters x t - - @variables c(..) - - ∂t = Differential(t) - ∂x = Differential(x) - ∂²x = Differential(x)^2 - - D₀ = 1.0 - R = 0.5 - cₑ = 2.5 - ℓ = 2.0 - α = 1 / 3 - - diff_eq = ∂t(c(x, t)) ~ ∂x((D₀ + α * c(x, t)) * ∂x(c(x, t))) - - bcs = [c(x, 0) ~ 0.0, # initial condition - R * (D₀ + α * c(0, t)) * ∂x(c(0, t)) ~ c(0, t) - cₑ, # Robin - ∂x(c(ℓ, t)) ~ 0.0] # no flux - - domains = [t ∈ Interval(0.0, 10.0), - x ∈ Interval(0.0, ℓ)] - - @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) - - # Method of lines discretization - dx = 0.1 - order = 2 - discretization = MOLFiniteDifference([x => dx], t, approx_order=order) - - # Convert the PDE problem into an ODE problem - prob = discretize(pdesys, discretization) - - sol = solve(prob, Rodas4()) -end - -@testset "Array u" begin - # Dependencies - N = 6 # number of dependent variables - - # Variables, parameters, and derivatives - @parameters x - @variables u[1:N](..) - Dx = Differential(x) - Dxx = Differential(x)^2 - - # Domain edges - x_min = 0.0 - x_max = 1.0 - - # Discretization parameters - dx = 0.1 - order = 2 - - #u = collect(u) - - # Equations - eqs = Vector{ModelingToolkit.Equation}(undef, N) - for i = 1:N - eqs[i] = Dxx(u[i](x)) ~ u[i](x) - end - - # Initial and boundary conditions - bcs = Vector{ModelingToolkit.Equation}(undef, 2 * N) - for i = 1:N - bcs[i] = Dx(u[i](x_min)) ~ 0.0 - end - - for i = 1:N - bcs[i+N] = u[i](x_max) ~ rand(StableRNG(0),) - end - - # Space and time domains - domains = [x ∈ Interval(x_min, x_max)] - - # PDE system - @named pdesys = PDESystem(eqs, bcs, domains, [x], collect([u[i](x) for i = 1:N])) - - # Method of lines discretization - discretization = MOLFiniteDifference([x => dx], nothing, approx_order=order) - prob = ModelingToolkit.discretize(pdesys, discretization) - - # # Solution of the ODE system - sol = NonlinearSolve.solve(prob, NewtonRaphson()) -end - -@testset "Non-Uniform Grid" begin - @parameters x y t - @variables u(..) v(..) - Dt = Differential(t) - Dx = Differential(x) - Dy = Differential(y) - Dxx = Differential(x)^2 - Dyy = Differential(y)^2 - - ∇²(u) = Dxx(u) + Dyy(u) - - brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.0 - - x_min = y_min = t_min = 0.0 - x_max = y_max = 1.0 - t_max = 11.5 - - α = 10.0 - - u0(x, y, t) = 22(y * (1 - y))^(3 / 2) - v0(x, y, t) = 27(x * (1 - x))^(3 / 2) - - eq = [Dt(u(x, y, t)) ~ 1.0 + Dy(v(x, y, t)) * u(x, y, t)^2 - 4.4 * u(x, y, t) + α * ∇²(u(x, y, t)) + brusselator_f(x, y, t), - Dt(v(x, y, t)) ~ -3.4 * Dx(u(x, y, t)) - v(x, y, t) * u(x, y, t)^2 + α * ∇²(v(x, y, t))] - - domains = [x ∈ Interval(x_min, x_max), - y ∈ Interval(y_min, y_max), - t ∈ Interval(t_min, t_max)] - - bcs = [u(x, y, 0) ~ u0(x, y, 0), - u(0, y, t) ~ 0.0, - u(1, y, t) ~ 0.0, - u(x, 0, t) ~ 0.0, - u(x, 1, t) ~ 0.0, - v(x, y, 0) ~ v0(x, y, 0), - v(0, y, t) ~ 0.0, - v(1, y, t) ~ 0.0, - v(x, 0, t) ~ 0.0, - v(x, 1, t) ~ 0.0] - - @named pdesys = PDESystem(eq, bcs, domains, [x, y, t], [u(x, y, t), v(x, y, t)]) - - # Method of lines discretization - N = 5 - - dx = [0.0, cumsum(rand(StableRNG(0), 0.05:0.01:0.2, 5))...] - if dx[end] < 1.0 - push!(dx, 1.0) - end - - dy = [0.0, cumsum(rand(StableRNG(0), 0.05:0.01:0.2, 5))...] - if dy[end] < 1.0 - push!(dy, 1.0) - end - - order = 2 - - discretization = MOLFiniteDifference([x => dx, y => dy], t, approx_order=order) - - prob = discretize(pdesys, discretization) - - sol = solve(prob, TRBDF2()) - -end - -@test_broken begin #@testset "2D variable connected to 1D variable at boundary #33" begin - @parameters t x r - @variables u(..) v(..) - Dt = Differential(t) - Dx = Differential(x) - Dxx = Differential(x)^2 - Dr = Differential(r) - Drr = Differential(r)^2 - - s = u(t, x) + v(t, x, 1) - - eqs = [Dt(u(t, x)) ~ Dxx(u(t, x)) + s, - Dt(v(t, x, r)) ~ Drr(v(t, x, r))] - bcs = [u(0, x) ~ 1, - v(0, x, r) ~ 1, - Dx(u(t, 0)) ~ 0.0, Dx(u(t, 1)) ~ 0.0, - Dr(v(t, x, 0)) ~ 0.0, Dr(v(t, x, 1)) ~ s] - - domains = [t ∈ Interval(0.0, 1.0), - x ∈ Interval(0.0, 1.0), - r ∈ Interval(0.0, 1.0)] - - @named pdesys = PDESystem(eqs, bcs, domains, [t, x, r], [u(t, x), v(t, x, r)]) - - # Method of lines discretization - dx = 0.1 - dr = 0.1 - order = 2 - discretization = MOLFiniteDifference([x => dx, r => dr], t, approx_order=order) - - # Convert the PDE problem into an ODE problem - prob = discretize(pdesys, discretization) - - sol = solve(prob, Tsit5()) -end diff --git a/test/pde_systems/MOLtest2_JSC.jl b/test/pde_systems/MOLtest2_JSC.jl deleted file mode 100644 index dadf55315..000000000 --- a/test/pde_systems/MOLtest2_JSC.jl +++ /dev/null @@ -1,478 +0,0 @@ -using ModelingToolkit, MethodOfLines, LinearAlgebra, OrdinaryDiffEq -using ModelingToolkit: operation, istree, arguments -using DomainSets -using NonlinearSolve -using StableRNGs -using JuliaSimCompiler -using Test - -# # Define some variables - - - -@testset "Testing discretization of varied systems" begin - @parameters x t - - @variables c(..) - - ∂t = Differential(t) - ∂x = Differential(x) - ∂²x = Differential(x)^2 - - D₀ = 1.5 - α = 0.15 - χ = 1.2 - R = 0.1 - cₑ = 2.0 - ℓ = 1.0 - Δx = 0.1 - - bcs = [ - # initial condition - c(x, 0) ~ 0.0, - # Robin BC - ∂x(c(0.0, t)) / (1 + exp(α * (c(0.0, t) - χ))) * R * D₀ + cₑ - c(0.0, t) ~ 0.0, - # no flux BC - ∂x(c(ℓ, t)) ~ 0.0] - - - # define space-time plane - domains = [x ∈ Interval(0.0, ℓ), t ∈ Interval(0.0, 5.0)] - - @test_broken begin #@testset "Test 01: ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))" begin - D = D₀ / (1.0 + exp(α * (c(x, t) - χ))) - diff_eq = ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t))) - @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) - discretization = MOLFiniteDifference([x => Δx], t) - prob = discretize(pdesys, discretization) - end - - @testset "Test 02: ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))" begin - D = 1.0 / (1.0 + exp(α * (c(x, t) - χ))) - diff_eq = ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t))) - @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) - discretization = MOLFiniteDifference([x => Δx], t) - prob = discretize(pdesys, discretization) - end - - @testset "Test 03: ∂t(c(x, t)) ~ ∂x(1.0 / (1.0/D₀ + exp(α * (c(x, t) - χ))/D₀) * ∂x(c(x, t)))" begin - diff_eq = ∂t(c(x, t)) ~ ∂x(1.0 / (1.0 / D₀ + exp(α * (c(x, t) - χ)) / D₀) * ∂x(c(x, t))) - @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) - discretization = MOLFiniteDifference([x => Δx], t) - prob = discretize(pdesys, discretization) - end - - @test_broken begin #@testset "Test 04: ∂t(c(x, t)) ~ ∂x(D₀ / (1.0 + exp(α * (c(x, t) - χ))) * ∂x(c(x, t)))" begin - diff_eq = ∂t(c(x, t)) ~ ∂x(D₀ / (1.0 + exp(α * (c(x, t) - χ))) * ∂x(c(x, t))) - @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) - discretization = MOLFiniteDifference([x => Δx], t) - prob = discretize(pdesys, discretization) - end - - @testset "Test 05: ∂t(c(x, t)) ~ ∂x(1/x * ∂x(c(x, t)))" begin - diff_eq = ∂t(c(x, t)) ~ ∂x(1 / x * ∂x(c(x, t))) - @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) - discretization = MOLFiniteDifference([x => Δx], t) - prob = discretize(pdesys, discretization) - end - - @test_broken begin #@testset "Test 06: ∂t(c(x, t)) ~ ∂x(x*∂x(c(x, t)))/c(x,t)" begin - diff_eq = ∂t(c(x, t)) ~ ∂x(x * ∂x(c(x, t))) / c(x, t) - @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) - discretization = MOLFiniteDifference([x => Δx], t) - prob = discretize(pdesys, discretization) - end - - @testset "Test 07: ∂t(c(x, t)) ~ ∂x(1/(1+c(x,t)) ∂x(c(x, t)))" begin - diff_eq = ∂t(c(x, t)) ~ ∂x(1 / (1 + c(x, t)) * ∂x(c(x, t))) - @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) - discretization = MOLFiniteDifference([x => Δx], t) - prob = discretize(pdesys, discretization) - end - - @test_broken begin #@testset "Test 08: ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x,t) * ∂x(c(x, t)))" begin - diff_eq = ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x, t) * ∂x(c(x, t))) - @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) - discretization = MOLFiniteDifference([x => Δx], t) - prob = discretize(pdesys, discretization) - end - - @test_broken begin #@testset "Test 09: ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x,t) * ∂x(c(x, t)))/(1+c(x,t))" begin - diff_eq = c(x, t) * ∂x(c(x, t) * ∂x(c(x, t))) / (1 + c(x, t)) ~ 0 - @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) - discretization = MOLFiniteDifference([x => Δx], t) - prob = discretize(pdesys, discretization) - end - - @test_broken begin #@testset "Test 10: ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x,t) * ∂x(c(x, t)))/(1+c(x,t))" begin - diff_eq = c(x, t) * ∂x(c(x, t)^(-1) * ∂x(c(x, t))) ~ 0 - @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) - discretization = MOLFiniteDifference([x => Δx], t) - prob = discretize(pdesys, discretization) - end - - @testset "Test 11: ∂t(c(x, t)) ~ ∂x(1/(1+c(x,t)^2) ∂x(c(x, t)))" begin - diff_eq = ∂t(c(x, t)) ~ ∂x(1 / (1 + c(x, t)^2) * ∂x(c(x, t))) - @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) - discretization = MOLFiniteDifference([x => Δx], t) - prob = discretize(pdesys, discretization) - end - -end - -@test_broken begin #@testset "Nonlinlap with flux interface boundary condition" begin - @parameters t x1 x2 - @variables c1(..) - @variables c2(..) - Dt = Differential(t) - - Dx1 = Differential(x1) - Dxx1 = Dx1^2 - - Dx2 = Differential(x2) - Dxx2 = Dx2^2 - - D1(c) = 1 + c / 10 - D2(c) = 1 / 10 + c / 10 - - eqs = [Dt(c1(t, x1)) ~ Dx1(D1(c1(t, x1)) * Dx1(c1(t, x1))), - Dt(c2(t, x2)) ~ Dx2(D2(c2(t, x2)) * Dx2(c2(t, x2)))] - - bcs = [c1(0, x1) ~ 1 + cospi(2 * x1), - c2(0, x2) ~ 1 + cospi(2 * x2), - Dx1(c1(t, 0)) ~ 0, - c1(t, 0.5) ~ c2(t, 0.5), - -D1(c1(t, 0.5)) * Dx1(c1(t, 0.5)) ~ -D2(c2(t, 0.5)) * Dx2(c2(t, 0.5)), - Dx2(c2(t, 1)) ~ 0] - - domains = [t ∈ Interval(0.0, 0.15), - x1 ∈ Interval(0.0, 0.5), - x2 ∈ Interval(0.5, 1.0)] - - @named pdesys = PDESystem(eqs, bcs, domains, - [t, x1, x2], [c1(t, x1), c2(t, x2)]) - - l = 40 - - disc = MOLFiniteDifference([x1 => l, x2 => l], t) - - prob = discretize(pdesys, disc) - - sol = solve(prob, FBDF(), saveat=0.01) - - x1_sol = sol[x1] - x2_sol = sol[x2] - t_sol = sol[t] - solc1 = sol[c1(t, x1)] - solc2 = sol[c2(t, x2)] - - solc = hcat(solc1[:, :], solc2[:, 2:end]) - - @test sol.retcode == SciMLBase.ReturnCode.Success -end - -@test_broken begin#@testset "Another boundaries appearing in equations case" begin - - g = 9.81 - - @parameters x z t - @variables φ(..) φ̃(..) η(..) - - Dt = Differential(t) - Dx = Differential(x) - Dz = Differential(z) - Dxx = Differential(x)^2 - Dzz = Differential(z)^2 - - eqs = [Dxx(φ(t, x, z)) + Dzz(φ(t, x, z)) ~ 0, - Dt(φ̃(t, x)) ~ -g * η(t, x), - Dt(η(t, x)) ~ Dz(φ(t, x, 1.0)) - ] - - bcs = [ - φ(0, x, z) ~ 0, - φ̃(0.0, x) ~ 0.0, - η(0.0, x) ~ cos(2 * π * x), - φ(t, x, 1.0) ~ φ̃(t, x), - Dx(φ(t, 0.0, z)) ~ 0.0, - Dx(φ(t, 1.0, z)) ~ 0.0, - Dz(φ(t, x, 0.0)) ~ 0.0, - Dx(φ̃(t, 0.0)) ~ 0.0, - Dx(φ̃(t, 1.0)) ~ 0.0, - Dx(η(t, 0.0)) ~ 0.0, - Dx(η(t, 1.0)) ~ 0.0, - ] - - domains = [x ∈ Interval(0.0, 1.0), - z ∈ Interval(0.0, 1.0), - t ∈ Interval(0.0, 10.0)] - - @named pdesys = PDESystem(eqs, bcs, domains, [t, x, z], - [φ(t, x, z), φ̃(t, x), η(t, x)]) - - - dx = 0.1 - dz = 0.1 - order = 2 - - discretization = MOLFiniteDifference([x => dx, z => dz], t, - approx_order=order, - grid_align=center_align) - - println("Discretization:") - prob = discretize(pdesys, discretization) -end - -@testset "Integrals in BCs" begin - β = 0.0005 - γ = 0.25 - amin = 0.0 - amax = 40.0 - - @parameters t a - @variables S(..) I(..) R(..) - Dt = Differential(t) - Da = Differential(a) - Ia = Integral(a in DomainSets.ClosedInterval(amin, amax)) - - - eqs = [Dt(S(t)) ~ -β * S(t) * Ia(I(a, t)), - Dt(I(a, t)) + Da(I(a, t)) ~ -γ * I(a, t), - Dt(R(t)) ~ γ * Ia(I(a, t))] - - bcs = [ - S(0) ~ 990.0, - I(0, t) ~ β * S(t) * Ia(I(a, t)), - I(a, 0) ~ 10.0 / 40.0, - R(0) ~ 0.0 - ] - - domains = [t ∈ (0.0, 40.0), a ∈ (0.0, 40.0)] - - @named pde_system = PDESystem(eqs, bcs, domains, [a, t], [S(t), I(a, t), R(t)]) - - da = 40 - discretization = MOLFiniteDifference([a => da], t) - - prob = MethodOfLines.discretize(pde_system, discretization) - - sol = solve(prob, FBDF()) - -end - -@test_broken begin #@testset "Dt in BCs" begin - # Parameters, variables, and derivatives - @parameters t x - @variables u(..) - Dt = Differential(t) - Dxx = Differential(x)^2 - - # 1D PDE and boundary conditions - eq = Dt(u(t, x)) ~ Dxx(u(t, x)) - bcs = [u(0, x) ~ 20, - Dt(u(t, 0)) ~ 100, # Heat source - Dt(u(t, 1)) ~ 0] # Zero flux - - # Space and time domains - domains = [t ∈ Interval(0.0, 1.0), - x ∈ Interval(0.0, 1.0)] - - # PDE system - @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)]) - - # Method of lines discretization - dx = 0.1 - order = 2 - discretization = MOLFiniteDifference([x => dx], t) - - # Convert the PDE problem into an ODE problem - prob = discretize(pdesys, discretization) - - # Solve ODE problem - sol = solve(prob, Rodas4(), saveat=0.2) - - discrete_x = sol[x] - discrete_t = sol[t] - solu = sol[u(t, x)] # Temperature should increase with time -end - -@test_broken begin #@testset "ODE connected to PDE at boundary" begin - @variables u(..) v(..) w(..) - @parameters t, r - Dt = Differential(t) - Dr = Differential(r) - Drr = Differential(r)^2 - - R = 1.0 - k₁ = 0.1 - k₂ = 0.1 - α = 1.0 - - u0 = 0.3 - v0 = 0.1 - w0 = 0.2 - - eq = [Dt(u(r, t)) ~ α * Drr(u(r, t)), - Dt(v(t)) ~ -k₁ * u(R, t) * v(t) + k₂ * w(t), - Dt(w(t)) ~ k₁ * u(R, t) * v(t) - k₂ * w(t) - ] - - bcs = [Dr(u(0, t)) ~ 0.0, - Dr(u(R, t)) ~ (-k₁ * u(R, t) * v(t) + k₂ * w(t)) / α, - u(r, 0) ~ u0, - v(0) ~ v0, - w(0) ~ w0 - ] - - domains = [t ∈ Interval(0.0, 10.0), - r ∈ Interval(0.0, R)] - - @named pdesys = PDESystem(eq, bcs, domains, [r, t], [u(r, t), v(t), w(t)]) - - dr = 0.1 - - disc = MOLFiniteDifference([r => dr], t) - - prob = discretize(pdesys, disc) - - sol = solve(prob, Rodas4P()) - - discrete_r = sol[r] - discrete_t = sol[t] - solu = sol[u(r, t)] - solv = sol[v(t)] - solw = sol[w(t)] -end - -@test_broken begin #@testset "ODE connected to PDE" begin - - @parameters t x - @variables u(..) v(..) - Dt = Differential(t) - Dx = Differential(x) - Dxx = Dx^2 - - # 1D PDE and boundary conditions - eqs = [Dt(u(t, x)) ~ Dxx(u(t, x)) + v(t), # This is the only line that is significantly changed from the test. - Dt(v(t)) ~ -v(t)] - - bcs = [u(0, x) ~ sin(x), - v(0) ~ 1, - u(t, 0) ~ 0, - Dx(u(t, 1)) ~ exp(-t) * cos(1)] - domains = [t ∈ Interval(0.0, 1.0), - x ∈ Interval(0.0, 1.0)] - @named pdesys = PDESystem(eqs, bcs, domains, [t, x], [u(t, x), v(t)]) - discretization = MOLFiniteDifference([x => 0.01], t) - prob = discretize(pdesys, discretization) - - sol = solve(prob, Tsit5()) - - discrete_x = sol[x] - discrete_t = sol[t] - solu = sol[u(t, x)] - solv = sol[v(t)] -end - -@test_broken begin #@testset "New style array variable conversion and interception" begin - # Parameters, variables, and derivatives - n_comp = 2 - @parameters t, x, p[1:n_comp], q[1:n_comp] - @variables u(..)[1:n_comp] - Dt = Differential(t) - Dx = Differential(x) - Dxx = Differential(x)^2 - params = Symbolics.scalarize(reduce(vcat,[p .=> [1.5, 2.0], q .=> [1.2, 1.8]])) - # 1D PDE and boundary conditions - - eqs = [Dt(u(t, x)[i]) ~ p[i] * Dxx(u(t, x)[i]) for i in 1:n_comp] - - bcs = [[u(0, x)[i] ~ q[i] * cos(x), - u(t, 0)[i] ~ sin(t), - u(t, 1)[i] ~ exp(-t) * cos(1), - Dx(u(t,0)[i]) ~ 0.0] for i in 1:n_comp] - bcs_collected = reduce(vcat, bcs) - - # Space and time domains - domains = [t ∈ Interval(0.0, 1.0), - x ∈ Interval(0.0, 1.0)] - - # PDE system - - @named pdesys = PDESystem(eqs, bcs_collected, domains, [t, x], [u(t, x)[i] for i in 1:n_comp], Symbolics.scalarize(params)) - - - # Method of lines discretization - dx = 0.1 - order = 2 - discretization = MOLFiniteDifference([x => dx], t; approx_order = order) - - # Convert the PDE problem into an ODE problem - prob = discretize(pdesys,discretization) #error occurs here - - # Solve ODE problem - sol = solve(prob, Tsit5(), saveat=0.2) - - # Test that the system is correctly constructed - varname1 = Symbol("u_Any[1]") - varname2 = Symbol("u_Any[2]") - - - vars = @variables $varname1(..), $varname2(..) - - @test sol[u(t, x)[1]] == sol[vars[1](t, x)] - @test sol[u(t, x)[2]] == sol[vars[2](t, x)] - - @test sol(0.1, 0.1, dv = u(t, x)[1]) == sol(0.1, 0.1, dv = vars[1](t, x)) - @test sol(0.1, 0.1, dv = u(t, x)[2]) == sol(0.1, 0.1, dv = vars[2](t, x)) -end - -@testset "Budkyo-Sellers, nonlinlap case" begin - using MethodOfLines, DomainSets, ModelingToolkit, OrdinaryDiffEq - - @parameters φ t - @variables Tₛ(..) - - Dt = Differential(t) - Dφ = Differential(φ) - - A = 210 # W/m^2 - B = 2 # W/m^2K - D = 0.6 # W/m^2K - - P₂(x) = 0.5*(3*(x)^2 - 1) - α_func(ϕ) = 0.354 + 0.25*P₂(sin(ϕ)) - Q_func(t, ϕ) =430*cos(ϕ) - - s_in_y = 31556952.0 - - φmin = -Float64(π/2)+0.1 - φmax = Float64(π/2)-0.1 - - # Water fraction - f = 0.7 - # rho is density of Water - ρ = 1025 # kg/m^3 - # specific heat of Water - c_w = 4186 #J/(kgK) - # Depth of water - H = 70 #m - - #Heat Cap of earth - C = f*ρ*c_w*H - - eq = Dt(Tₛ(t, φ)) ~ ((1 - α_func(φ))*Q_func(t, φ) - A - B*Tₛ(t, φ) + D*Dφ(cos(φ)*Dφ(Tₛ(t, φ)))/cos(φ))/C - - bcs = [Tₛ(0, φ) ~ 12 - 40*P₂(sin(φ)), - Dφ(Tₛ(t, φmin)) ~ 0.0, - Dφ(Tₛ(t, φmax)) ~ 0.0] - domains = [t in IntervalDomain(0, 19*s_in_y), φ in IntervalDomain(φmin, φmax)] - - @named sys = PDESystem(eq, bcs, domains, [t, φ], [Tₛ(t, φ)]) - - discretization = MOLFiniteDifference([φ => 80], t, order = 4) - - prob = discretize(sys, discretization) # ERROR HERE - - sol = solve(prob, FBDF(), saveat = s_in_y) -end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 7c3e8be40..39e3b6b89 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -128,17 +128,4 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") include("pde_systems/MOL_1D_Linear_Convection.jl") end end - ############### JSC ##################### - if GROUP == "All" || GROUP == "MOL_Interface2_JSC" - Pkg.add("JuliaSimCompiler") - @time @safetestset "MOLFiniteDifference Interface" begin - include("pde_systems/MOLtest2_JSC.jl") - end - end - if GROUP == "All" || GROUP == "MOL_Interface1_JSC" - Pkg.add("JuliaSimCompiler") - @time @safetestset "MOLFiniteDifference Interface" begin - include("pde_systems/MOLtest1_JSC.jl") - end - end end From 3f71d8244e3979aface462707c851ee66b7aa850 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Wed, 20 Sep 2023 20:42:16 +0100 Subject: [PATCH 33/48] Update Project.toml --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index affa1b02b..fcae0aa16 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MethodOfLines" uuid = "94925ecb-adb7-4558-8ed8-f975c56a0bf4" -authors = ["Alex Jones, "] -version = "0.10.1" +authors = ["Alex Jones, "] +version = "0.10.0" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" From 3139607dc0e4fdccf5b63e5504c85e2b3e2dcd48 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Wed, 27 Sep 2023 16:02:33 +0100 Subject: [PATCH 34/48] Update Project.toml --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index fcae0aa16..96397a032 100644 --- a/Project.toml +++ b/Project.toml @@ -33,6 +33,7 @@ ModelingToolkit = "8" OrdinaryDiffEq = "6" PDEBase = "0.1.4" PrecompileTools = "1" +RuntimeGeneratedFunctions = "0.5.9" SafeTestsets = "0.0.1" SciMLBase = "1.94, 2" StaticArrays = "1" From 93ed1b8fafe47c99b3f39d4b25f4df81bb1ffbcc Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Wed, 27 Sep 2023 15:34:28 +0000 Subject: [PATCH 35/48] CompatHelper: bump compat for MethodOfLines to 0.10 for package docs, (keep existing compat) --- docs/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index e1642e490..e208b9d75 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -12,7 +12,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" DifferentialEquations = "7.6" Documenter = "1" DomainSets = "0.5, 0.6" -MethodOfLines = "0.9" +MethodOfLines = "0.9, 0.10" ModelingToolkit = "8" NonlinearSolve = "1" OrdinaryDiffEq = "6" From 27c1bd27fed028cb781d9e42b7fa0fcecc4f1dd7 Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Fri, 29 Sep 2023 18:07:54 +0000 Subject: [PATCH 36/48] CompatHelper: bump compat for NonlinearSolve to 2 for package docs, (keep existing compat) --- docs/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index e1642e490..6daaac5ac 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -14,6 +14,6 @@ Documenter = "1" DomainSets = "0.5, 0.6" MethodOfLines = "0.9" ModelingToolkit = "8" -NonlinearSolve = "1" +NonlinearSolve = "1, 2" OrdinaryDiffEq = "6" Plots = "1" From 67a5b3d43a8709d3d61521edcf0f84629d2b196a Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Fri, 6 Oct 2023 06:43:08 +0000 Subject: [PATCH 37/48] docs: fix rendering of the plots in the heat equation example --- docs/src/tutorials/heat.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/src/tutorials/heat.md b/docs/src/tutorials/heat.md index d817351d8..aa6eb73b8 100644 --- a/docs/src/tutorials/heat.md +++ b/docs/src/tutorials/heat.md @@ -48,11 +48,11 @@ solu = sol[u(t, x)] using Plots plt = plot() -for i in 1:length(discrete_t) +for i in eachindex(discrete_t) plot!(discrete_x, solu[i, :], label="Numerical, t=$(discrete_t[i])") scatter!(discrete_x, u_exact(discrete_x, discrete_t[i]), label="Exact, t=$(discrete_t[i])") end -display(plt) +plt ``` ### Neumann boundary conditions @@ -103,11 +103,11 @@ solu = sol[u(t, x)] using Plots plt = plot() -for i in 1:length(discrete_t) +for i in eachindex(discrete_t) plot!(discrete_x, solu[i, :], label="Numerical, t=$(discrete_t[i])") scatter!(discrete_x, u_exact(discrete_x, discrete_t[i]), label="Exact, t=$(discrete_t[i])") end -display(plt) +plt ``` ### Robin boundary conditions @@ -159,11 +159,11 @@ solu = sol[u(t, x)] using Plots plt = plot() -for i in 1:length(discrete_t) +for i in eachindex(discrete_t) plot!(discrete_x, solu[i, :], label="Numerical, t=$(discrete_t[i])") scatter!(discrete_x, u_exact(discrete_x, discrete_t[i]), label="Exact, t=$(discrete_t[i])") end -display(plt) +plt ``` From e3a056f1a8f4db1a147bdb47081df2a4f0a729be Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Wed, 11 Oct 2023 11:58:36 +0100 Subject: [PATCH 38/48] fix fail --- test/runtests.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 7acef7ba8..257e158fd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -30,7 +30,6 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") @time @safetestset "MOLFiniteDifference Interface: Advanced Nonlinear Diffusion" begin include("pde_systems/nonlinear_laplacian_advanced.jl") end - end end if GROUP == "All" || GROUP == "Sol_Interface" From 261972722071f63fb778ca492337c86b37071e66 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Wed, 11 Oct 2023 12:06:04 +0100 Subject: [PATCH 39/48] fix typo --- test/pde_systems/MOL_1D_Linear_Diffusion.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/pde_systems/MOL_1D_Linear_Diffusion.jl b/test/pde_systems/MOL_1D_Linear_Diffusion.jl index f4ff722b4..09f01484b 100644 --- a/test/pde_systems/MOL_1D_Linear_Diffusion.jl +++ b/test/pde_systems/MOL_1D_Linear_Diffusion.jl @@ -42,7 +42,7 @@ using ModelingToolkit: Differential for disc in [discretization, discretization_edge, discretization_centered, discretization_approx_order4] # Convert the PDE problem into an ODE problem - prob = discrebranch tize(pdesys, disc) + prob = discretize(pdesys, disc) # Solve ODE problem # Solve ODE problem sol = solve(prob, Tsit5(), saveat=0.1) From 41ad19426eed0ac14d9940f3cb62f1d88a5ced1f Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Wed, 11 Oct 2023 12:46:19 +0100 Subject: [PATCH 40/48] fix typo --- .../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 b8e99144b..8325c8794 100644 --- a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl +++ b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl @@ -116,7 +116,7 @@ end 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, indexmap, bcmap, 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, indexmapp, bcmap, 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, indexmap, bcmap, 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))(*(~~a, $(Differential(x))(u), ~~d))), ~~c), ~e) => /(*(~b..., ~c..., cartesian_nonlinear_laplacian(*(~a..., ~d...), Idx(II, s, u, indexmap), derivweights, s, indexmap, bcmap, depvars, x, u)), replacevals(~e, s, u, depvars, II, indexmap)) for x in ivs(u, s)]) for u in depvars], init = [])) From e2974eaffa5d7810e87859d4bc99ca8a966b3be7 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Wed, 11 Oct 2023 14:15:50 +0100 Subject: [PATCH 41/48] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 96397a032..e897ccd0f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MethodOfLines" uuid = "94925ecb-adb7-4558-8ed8-f975c56a0bf4" authors = ["Alex Jones, "] -version = "0.10.0" +version = "0.10.1" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" From f4139aabcd1e8ba2c3a047906f2700bad41f0f1b Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Wed, 18 Oct 2023 15:40:01 +0100 Subject: [PATCH 42/48] a --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index c4b4b9719..2db1a2f8d 100644 --- a/.gitignore +++ b/.gitignore @@ -16,3 +16,4 @@ dev/* # Build artifacts for creating documentation generated by the Documenter package docs/build/ docs/site/ +.DS_Store \ No newline at end of file From 422e0a0daf3489cc092946e8c539b0d8b6df9a68 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Thu, 19 Oct 2023 17:21:57 +0100 Subject: [PATCH 43/48] add callbacks --- src/MethodOfLines.jl | 4 ++- .../differential_discretizer.jl | 4 ++- .../generate_finite_difference_rules.jl | 4 ++- .../schemes/callbacks/callback_rules.jl | 5 +++ src/interface/MOLFiniteDifference.jl | 5 +-- src/interface/callbacks.jl | 13 +++++++ test/components/callbacks.jl | 35 +++++++++++++++++++ 7 files changed, 65 insertions(+), 5 deletions(-) create mode 100644 src/discretization/schemes/callbacks/callback_rules.jl create mode 100644 src/interface/callbacks.jl create mode 100644 test/components/callbacks.jl diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl index 7a818230e..5bd8e6107 100644 --- a/src/MethodOfLines.jl +++ b/src/MethodOfLines.jl @@ -50,6 +50,7 @@ import Base.ndims # Interface include("interface/grid_types.jl") include("interface/scheme_types.jl") +include("interface/callbacks.jl") include("interface/disc_strategy_types.jl") include("interface/MOLFiniteDifference.jl") @@ -72,6 +73,7 @@ include("discretization/schemes/upwind_difference/upwind_diff_weights.jl") include("discretization/schemes/half_offset_weights.jl") include("discretization/schemes/extrapolation_weights.jl") include("discretization/differential_discretizer.jl") +include("discretization/schemes/callbacks/callback_rules.jl") # System Parsing include("system_parsing/pde_system_transformation.jl") @@ -104,6 +106,6 @@ include("precompile.jl") # Export export MOLFiniteDifference, discretize, symbolic_discretize, ODEFunctionExpr, generate_code, grid_align, edge_align, center_align, get_discrete, chebyspace -export UpwindScheme, WENOScheme, FunctionalScheme +export UpwindScheme, WENOScheme, FunctionalScheme, MOLDiscCallback end diff --git a/src/discretization/differential_discretizer.jl b/src/discretization/differential_discretizer.jl index 5a6250c76..d23615d0e 100644 --- a/src/discretization/differential_discretizer.jl +++ b/src/discretization/differential_discretizer.jl @@ -8,6 +8,7 @@ struct DifferentialDiscretizer{T, D1, S} <: AbstractDifferentialDiscretizer interpmap::Dict{Num,DerivativeOperator} orders::Dict{Num,Vector{Int}} boundary::Dict{Num,DerivativeOperator} + callbacks end function PDEBase.construct_differential_discretizer(pdesys, s::DiscreteSpace, discretization::MOLFiniteDifference, orders) @@ -15,6 +16,7 @@ function PDEBase.construct_differential_discretizer(pdesys, s::DiscreteSpace, di bcs = pdesys.bcs isa Vector ? pdesys.bcs : [pdesys.bcs] approx_order = discretization.approx_order advection_scheme = discretization.advection_scheme + callbacks = discretization.callbacks upwind_order = advection_scheme isa UpwindScheme ? advection_scheme.order : 1 differentialmap = Array{Pair{Num,DerivativeOperator},1}() @@ -63,5 +65,5 @@ function PDEBase.construct_differential_discretizer(pdesys, s::DiscreteSpace, di push!(boundary, x => BoundaryInterpolatorExtrapolator(max(6, approx_order), dx)) end - return DifferentialDiscretizer{eltype(orders),typeof(Dict(differentialmap)),typeof(advection_scheme)}(approx_order, advection_scheme, Dict(differentialmap), (Dict(nonlinlap_inner), Dict(nonlinlap_outer)), (Dict(windpos), Dict(windneg)), Dict(interp), Dict(orders), Dict(boundary)) + return DifferentialDiscretizer{eltype(orders),typeof(Dict(differentialmap)),typeof(advection_scheme)}(approx_order, advection_scheme, Dict(differentialmap), (Dict(nonlinlap_inner), Dict(nonlinlap_outer)), (Dict(windpos), Dict(windneg)), Dict(interp), Dict(orders), Dict(boundary), callbacks) end diff --git a/src/discretization/generate_finite_difference_rules.jl b/src/discretization/generate_finite_difference_rules.jl index 8d2ad1b5d..e5eacab33 100644 --- a/src/discretization/generate_finite_difference_rules.jl +++ b/src/discretization/generate_finite_difference_rules.jl @@ -80,6 +80,8 @@ function generate_finite_difference_rules(II::CartesianIndex, s::DiscreteSpace, integration_rules = [] end + cb_rules = generate_cb_rules(II, s, depvars, derivweights, bmap, indexmap, terms) + integration_rules = vcat(integration_rules, vec(generate_whole_domain_integration_rules(II, s, depvars, indexmap, terms))) - return vcat(vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(mixed_deriv_rules_cartesian), vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules) + return vcat(vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(mixed_deriv_rules_cartesian), vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules, cb_rules) end diff --git a/src/discretization/schemes/callbacks/callback_rules.jl b/src/discretization/schemes/callbacks/callback_rules.jl new file mode 100644 index 000000000..cc88217e5 --- /dev/null +++ b/src/discretization/schemes/callbacks/callback_rules.jl @@ -0,0 +1,5 @@ +function generate_cb_rules(II, s, depvars, derivweights, bmap, indexmap, terms) + cbs = derivweights.callbacks + cb_rules = [cb.sym => cb(s) for cb in cbs] + return cb_rules +end \ No newline at end of file diff --git a/src/interface/MOLFiniteDifference.jl b/src/interface/MOLFiniteDifference.jl index a4a7d7173..968e105e2 100644 --- a/src/interface/MOLFiniteDifference.jl +++ b/src/interface/MOLFiniteDifference.jl @@ -33,11 +33,12 @@ struct MOLFiniteDifference{G,D} <: AbstractEquationSystemDiscretization use_ODAE::Bool disc_strategy::D useIR::Bool + callbacks kwargs end # Constructors. If no order is specified, both upwind and centered differences will be 2nd order -function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, advection_scheme = UpwindScheme(), grid_align=CenterAlignedGrid(), discretization_strategy = ScalarizedDiscretization(), upwind_order = nothing, should_transform = true, use_ODAE = false, useIR = true, kwargs...) +function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, advection_scheme = UpwindScheme(), grid_align=CenterAlignedGrid(), discretization_strategy = ScalarizedDiscretization(), upwind_order = nothing, should_transform = true, use_ODAE = false, useIR = true, callbacks = [], kwargs...) if upwind_order !== nothing @warn "`upwind_order` no longer does anything, and will be removed in a future release. See the docs for the current interface." end @@ -50,7 +51,7 @@ function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, advection_sche dxs = dxs isa Dict ? dxs : Dict(dxs) - return MOLFiniteDifference{typeof(grid_align), typeof(discretization_strategy)}(dxs, time, approx_order, advection_scheme, grid_align, should_transform, use_ODAE, discretization_strategy, useIR, kwargs) + return MOLFiniteDifference{typeof(grid_align), typeof(discretization_strategy)}(dxs, time, approx_order, advection_scheme, grid_align, should_transform, use_ODAE, discretization_strategy, useIR, callbacks, kwargs) end PDEBase.get_time(disc::MOLFiniteDifference) = disc.time diff --git a/src/interface/callbacks.jl b/src/interface/callbacks.jl new file mode 100644 index 000000000..2683e1a4f --- /dev/null +++ b/src/interface/callbacks.jl @@ -0,0 +1,13 @@ +struct MOLDiscCallback + f::Function + p + sym +end + +function MOLDiscCallback(s, disc_ps) + symh = hash(f) + sym = unwrap(first(@parameters(Symbol("MOLDiscCallback_$symh")))) + return MOLDiscCallback(f, disc_ps, sym) +end + +(M::MOLDiscCallback)(s) = M.f(s, M.p) \ No newline at end of file diff --git a/test/components/callbacks.jl b/test/components/callbacks.jl new file mode 100644 index 000000000..1a8721217 --- /dev/null +++ b/test/components/callbacks.jl @@ -0,0 +1,35 @@ +using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils + +@testset "Discrete callback" begin + @parameters x, t + @variables u(..) + Dt = Differential(t) + Dxx = Differential(x)^2 + t_min = 0.0 + t_max = 2.0 + x_min = 0.0 + x_max = 2.0 + + a = 0.1*pi + + eq = Dt(u(t, x)) ~ a*Dxx(u(t, x)) + + bcs = [u(t_min, x) ~ sinpi(x), u(t, x_min) ~ sinpi(t+x), u(t, x_max) ~ sinpi(t+x)] + + domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)] + + pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)]) + + cb = MOLDiscCallback((s, p) -> prod(p), [0.1, pi]) + a_cb = cb.sym + cbeq = Dt(u(t, x)) ~ a_cb*Dxx(u(t, x)) + + cbpdesys = PDESystem(cbeq, bcs, domains, [t, x], [u(t, x)], callbacks = [cb]) + + disc = MOLFiniteDifference([x => 0.1], t; approx_order=2) + + prob = discretize(pdesys, disc) + cbprob = discretize(cbpdesys, disc) + + @test solve(prob, Tsit5(), saveat=0.1) ≈ solve(cbprob, Tsit5(), saveat=0.1) +end \ No newline at end of file From ac10b6af34e78726fc9577fd0d8f4bc77ff19c20 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Thu, 19 Oct 2023 17:38:43 +0100 Subject: [PATCH 44/48] runtests --- src/discretization/generate_finite_difference_rules.jl | 2 +- test/runtests.jl | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/discretization/generate_finite_difference_rules.jl b/src/discretization/generate_finite_difference_rules.jl index e5eacab33..ac11911ec 100644 --- a/src/discretization/generate_finite_difference_rules.jl +++ b/src/discretization/generate_finite_difference_rules.jl @@ -83,5 +83,5 @@ function generate_finite_difference_rules(II::CartesianIndex, s::DiscreteSpace, cb_rules = generate_cb_rules(II, s, depvars, derivweights, bmap, indexmap, terms) integration_rules = vcat(integration_rules, vec(generate_whole_domain_integration_rules(II, s, depvars, indexmap, terms))) - return vcat(vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(mixed_deriv_rules_cartesian), vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules, cb_rules) + return vcat(cb_rules, vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(mixed_deriv_rules_cartesian), vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules) end diff --git a/test/runtests.jl b/test/runtests.jl index 257e158fd..263692a8f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -71,6 +71,9 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") @time @safetestset "ODEFunction" begin include("components/ODEFunction_test.jl") end + @time @safetestset "Discrete Callbacks" begin + include("components/callbacks.jl") + end #@time @safetestset "Finite Difference Schemes" begin include("components/finite_diff_schemes.jl") end end From ccc47f9fd53aec7c3b02b10ccd38818da239f6ba Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Tue, 14 Nov 2023 18:47:46 +0000 Subject: [PATCH 45/48] callbacks passing --- src/interface/callbacks.jl | 2 +- src/interface/solution/timeindep.jl | 2 +- src/system_parsing/pde_system_transformation.jl | 3 +++ test/components/callbacks.jl | 13 ++++++++----- 4 files changed, 13 insertions(+), 7 deletions(-) diff --git a/src/interface/callbacks.jl b/src/interface/callbacks.jl index 2683e1a4f..5c5fc427e 100644 --- a/src/interface/callbacks.jl +++ b/src/interface/callbacks.jl @@ -4,7 +4,7 @@ struct MOLDiscCallback sym end -function MOLDiscCallback(s, disc_ps) +function MOLDiscCallback(f, disc_ps) symh = hash(f) sym = unwrap(first(@parameters(Symbol("MOLDiscCallback_$symh")))) return MOLDiscCallback(f, disc_ps, sym) diff --git a/src/interface/solution/timeindep.jl b/src/interface/solution/timeindep.jl index 8e3374952..b1239c845 100644 --- a/src/interface/solution/timeindep.jl +++ b/src/interface/solution/timeindep.jl @@ -11,7 +11,7 @@ function SciMLBase.PDENoTimeSolution(sol::SciMLBase.NonlinearSolution{T}, metada umap = mapreduce(vcat, dvs) do u let discu = discretespace.discvars[u] solu = map(CartesianIndices(discu)) do I - i = sym_to_index(discu[I], odesys.states) + i = sym_to_index(discu[I], get_states(odesys)) # Handle Observed if i !== nothing sol.u[i] diff --git a/src/system_parsing/pde_system_transformation.jl b/src/system_parsing/pde_system_transformation.jl index 681b00866..ca722b151 100644 --- a/src/system_parsing/pde_system_transformation.jl +++ b/src/system_parsing/pde_system_transformation.jl @@ -35,6 +35,9 @@ function PDEBase.transform_pde_system!(v::PDEBase.VariableMap, boundarymap, sys: end function PDEBase.should_transform(pdesys::PDESystem, disc::MOLFiniteDifference, boundarymap) + if !disc.should_transform + return false + end if has_interfaces(boundarymap) @warn "The system contains interface boundaries, which are not compatible with system transformation. The system will not be transformed. Please post an issue if you need this feature." return false diff --git a/test/components/callbacks.jl b/test/components/callbacks.jl index 1a8721217..a1d1236f9 100644 --- a/test/components/callbacks.jl +++ b/test/components/callbacks.jl @@ -1,4 +1,4 @@ -using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils +using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils, OrdinaryDiffEq @testset "Discrete callback" begin @parameters x, t @@ -18,18 +18,21 @@ using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)] - pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)]) + @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)]) cb = MOLDiscCallback((s, p) -> prod(p), [0.1, pi]) a_cb = cb.sym cbeq = Dt(u(t, x)) ~ a_cb*Dxx(u(t, x)) - cbpdesys = PDESystem(cbeq, bcs, domains, [t, x], [u(t, x)], callbacks = [cb]) + @named cbpdesys = PDESystem(cbeq, bcs, domains, [t, x], [u(t, x)]) - disc = MOLFiniteDifference([x => 0.1], t; approx_order=2) + disc = MOLFiniteDifference([x => 0.1], t; approx_order=2, callbacks = [cb]) prob = discretize(pdesys, disc) cbprob = discretize(cbpdesys, disc) - @test solve(prob, Tsit5(), saveat=0.1) ≈ solve(cbprob, Tsit5(), saveat=0.1) + sol1 = solve(prob, Tsit5(), saveat=0.1) + sol2 = solve(cbprob, Tsit5(), saveat=0.1) + + @test sol1[u(t, x)] ≈ sol2[u(t, x)] end \ No newline at end of file From 135c86a51f6b2a75aa01f9129d81c8bfa94c0c83 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Tue, 14 Nov 2023 19:55:15 +0000 Subject: [PATCH 46/48] test --- test/components/callbacks.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/components/callbacks.jl b/test/components/callbacks.jl index a1d1236f9..e5a6b89ea 100644 --- a/test/components/callbacks.jl +++ b/test/components/callbacks.jl @@ -5,7 +5,7 @@ using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils @variables u(..) Dt = Differential(t) Dxx = Differential(x)^2 - t_min = 0.0 + t_min = 0.0 t_max = 2.0 x_min = 0.0 x_max = 2.0 From 7966aa2e15b33730135e6ed2a32c26da64ff852e Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Tue, 14 Nov 2023 19:55:59 +0000 Subject: [PATCH 47/48] update compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e897ccd0f..c1b21a2a8 100644 --- a/Project.toml +++ b/Project.toml @@ -31,7 +31,7 @@ Interpolations = "0.14" Latexify = "0.15, 0.16" ModelingToolkit = "8" OrdinaryDiffEq = "6" -PDEBase = "0.1.4" +PDEBase = "0.1.6" PrecompileTools = "1" RuntimeGeneratedFunctions = "0.5.9" SafeTestsets = "0.0.1" From a04e094fd64cc8dba8b3b32fc3bcb26485b37234 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Tue, 14 Nov 2023 22:10:34 +0000 Subject: [PATCH 48/48] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index c1b21a2a8..464e8a9cf 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MethodOfLines" uuid = "94925ecb-adb7-4558-8ed8-f975c56a0bf4" authors = ["Alex Jones, "] -version = "0.10.1" +version = "0.10.2" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"