diff --git a/src/discretization/generate_finite_difference_rules.jl b/src/discretization/generate_finite_difference_rules.jl index 7ae8b4760..f4ce35b61 100644 --- a/src/discretization/generate_finite_difference_rules.jl +++ b/src/discretization/generate_finite_difference_rules.jl @@ -77,24 +77,29 @@ function generate_finite_difference_rules( # Spherical diffusion scheme spherical_diffusion_rules = generate_spherical_diffusion_rules( II, s, depvars, derivweights, bmap, indexmap, split_additive_terms(pde)) - integration_rules = vec(generate_euler_integration_rules( - II, s, depvars, indexmap, terms)) + + # Integration schemes for boundaries: [lower-bound, x] and [x, upper-bound] + integration_rules_from_boundary_to_x = generate_euler_integration_rules( + II, s, depvars, indexmap, terms) else central_deriv_rules_cartesian = [] advection_rules = [] nonlinlap_rules = [] spherical_diffusion_rules = [] mixed_deriv_rules_cartesian = [] - integration_rules = [] + integration_rules_from_boundary_to_x = [] 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))) + # Integration schemes for boundaries: [lower-bound, upper-bound] + whole_domain_integration_rules = generate_whole_domain_integration_rules( + II, s, depvars, indexmap, terms) + 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) + vec(advection_rules), + vec(integration_rules_from_boundary_to_x), vec(whole_domain_integration_rules)) end function generate_finite_difference_rules( @@ -112,10 +117,12 @@ function generate_finite_difference_rules( advection_rules = [] nonlinlap_rules = [] spherical_diffusion_rules = [] - integration_rules = [] + integration_rules_from_boundary_to_x = [] + + whole_domain_integration_rules = generate_whole_domain_integration_rules( + II, s, depvars, 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(central_deriv_rules_cartesian), vec(advection_rules), integration_rules) + vec(central_deriv_rules_cartesian), vec(advection_rules), + vec(integration_rules_from_boundary_to_x), vec(whole_domain_integration_rules)) end diff --git a/src/discretization/schemes/integral_expansion/integral_expansion.jl b/src/discretization/schemes/integral_expansion/integral_expansion.jl index 1b8e777ed..3949f83e6 100644 --- a/src/discretization/schemes/integral_expansion/integral_expansion.jl +++ b/src/discretization/schemes/integral_expansion/integral_expansion.jl @@ -1,66 +1,90 @@ # use the trapezoid rule -function _euler_integral(II, s, jx, u, ufunc, dx::Number) #where {T,N,Wind,DX<:Number} +function _euler_integral_lower_to_II(II, s, jx, u, ufunc, dx) #where {T,N,Wind,DX<:Number} j, x = jx - if II[j] == 1 + if II[j] == 1 # recursively arrived at lower end of the domain return Num(0) end # unit index in direction of the derivative I1 = unitindex(ndims(u, s), j) # dx for multiplication Itap = [II - I1, II] - weights = [dx / 2, dx / 2] + if dx isa Number # Uniform dx + weights = [dx / 2, dx / 2] + elseif dx isa AbstractVector # Nonuniform dx + weights = fill(dx[II[j] - 1] / 2, 2) + else + error("Unsupported type of dx: $(type(dx))") + end + + # sym_do computes from II to II - I1, + # and recursive call computes from II - I1 to lower end of domain return sym_dot(weights, ufunc(u, Itap, x)) + - _euler_integral(II - I1, s, jx, u, ufunc, dx) + _euler_integral_lower_to_II(II - I1, s, jx, u, ufunc, dx) end - -# Nonuniform dx -function _euler_integral(II, s, jx, u, ufunc, dx::AbstractVector) #where {T,N,Wind,DX<:Number} +function _euler_integral_II_to_upper(II, s, jx, u, ufunc, dx) #where {T,N,Wind,DX<:Number} j, x = jx - if II[j] == 1 + if II[j] == length(s, x) # recursively arrived at upper end of the domain return Num(0) end # unit index in direction of the derivative I1 = unitindex(ndims(u, s), j) # dx for multiplication - Itap = [II - I1, II] - weights = fill(dx[II[j] - 1] / 2, 2) + Itap = [II, II + I1] + if dx isa Number # Uniform dx + weights = [dx / 2, dx / 2] + elseif dx isa AbstractVector # Nonuniform dx + weights = fill(dx[II[j] + 1] / 2, 2) + else + error("Unsupported type of dx: $(type(dx))") + end + + # sym_do computes from II to II + I1, + # and recursive call computes from II + I1 to upper end of domain return sym_dot(weights, ufunc(u, Itap, x)) + - _euler_integral(II - I1, s, jx, u, ufunc, dx) + _euler_integral_II_to_upper(II + I1, s, jx, u, ufunc, dx) end -function euler_integral(II, s, jx, u, ufunc) +# An integral from II to end of domain +function euler_integral(method::Symbol, II, s, jx, u, ufunc) j, x = jx dx = s.dxs[x] - return _euler_integral(II, s, jx, u, ufunc, dx) + if method == :lower_boundary_to_x + return _euler_integral_lower_to_II(II, s, jx, u, ufunc, dx) + elseif method == :x_to_upper_boundary + return _euler_integral_II_to_upper(II, s, jx, u, ufunc, dx) + # error("Method :x_to_upper_boundary is not implemented for euler_integral.") + else + error("Unsupported method for euler_integral: $method") + end end # An integral across the whole domain (xmin .. xmax) function whole_domain_integral(II, s, jx, u, ufunc) j, x = jx dx = s.dxs[x] - if II[j] == length(s, x) - return _euler_integral(II, s, jx, u, ufunc, dx) - end - - dist2max = length(s, x) - II[j] + dist2max = length(s, x) - II[j] # distance to the end of the domain I1 = unitindex(ndims(u, s), j) Imax = II + dist2max * I1 - return _euler_integral(Imax, s, jx, u, ufunc, dx) + + return _euler_integral_lower_to_II(Imax, s, jx, u, ufunc, dx) end @inline function generate_euler_integration_rules( II::CartesianIndex, s::DiscreteSpace, depvars, indexmap, terms) ufunc(u, I, x) = s.discvars[u][I] - eulerrules = reduce(safe_vcat, - [[Integral(x in DomainSets.ClosedInterval(s.vars.intervals[x][1], - Num(x)))(u) => euler_integral( - Idx(II, s, u, indexmap), s, (x2i(s, u, x), x), u, ufunc) - for x in ivs(u, s)] - for u in depvars], - init = []) + eulerrules = reduce(safe_vcat, [ + reduce(safe_vcat, [ + [ + # integrals from lower domain end to x: + Integral(x in DomainSets.ClosedInterval(s.vars.intervals[x][1], Num(x)))(u) => euler_integral(:lower_boundary_to_x, Idx(II, s, u, indexmap), s, (x2i(s, u, x), x), u, ufunc), + # integrals from x to upper domain end: + Integral(x in DomainSets.ClosedInterval(Num(x), s.vars.intervals[x][2]))(u) => euler_integral(:x_to_upper_boundary, Idx(II, s, u, indexmap), s, (x2i(s, u, x), x), u, ufunc) + ] + for x in ivs(u, s)]) for u in depvars] ) + return eulerrules end diff --git a/src/system_parsing/pde_system_transformation.jl b/src/system_parsing/pde_system_transformation.jl old mode 100644 new mode 100755 index 9c47960e5..78dbde0df --- a/src/system_parsing/pde_system_transformation.jl +++ b/src/system_parsing/pde_system_transformation.jl @@ -156,11 +156,13 @@ function descend_to_incompatible(term, v) end elseif op isa Integral if any(isequal(op.domain.variables), v.x̄) - euler = isequal( - op.domain.domain.left, v.intervals[op.domain.variables][1]) && - isequal(op.domain.domain.right, Num(op.domain.variables)) - whole = isequal( - op.domain.domain.left, v.intervals[op.domain.variables][1]) && + euler = ( # from lower domain boundary to x: + isequal(op.domain.domain.left, v.intervals[op.domain.variables][1]) && + isequal(op.domain.domain.right, Num(op.domain.variables))) || + (# from x to upper domain boundary: + isequal(op.domain.domain.left, Num(op.domain.variables)) && + isequal(op.domain.domain.right, v.intervals[op.domain.variables][2])) + whole = isequal(op.domain.domain.left, v.intervals[op.domain.variables][1]) && isequal(op.domain.domain.right, v.intervals[op.domain.variables][2]) if any([euler, whole]) u = arguments(term)[1] @@ -168,10 +170,10 @@ function descend_to_incompatible(term, v) @assert out==(nothing, false) "Integral $term must be purely of a variable, got $u. Try wrapping the integral argument with an auxiliary variable." return (nothing, nothing, false) else - throw(ArgumentError("Integration Domain only supported for integrals from start of iterval to the variable, got $(op.domain.domain) in $(term)")) + throw(ArgumentError("Integration Domain only supported for integrals across whole interval, or from an interval boundary (upper or lower) to the independent variable, got $(op.domain.domain) in $(term)")) end else - throw(ArgumentError("Integral must be with respect to the independent variable in its upper-bound, got $(op.domain.variables) in $(term)")) + throw(ArgumentError("Integral must be with respect to the independent variable in its upper- or lower-bound, got $(op.domain.variables) in $(term)")) end end for arg in SU.arguments(term) diff --git a/test/pde_systems/MOL_1D_Integration.jl b/test/pde_systems/MOL_1D_Integration.jl index aef8bab5d..b66217780 100644 --- a/test/pde_systems/MOL_1D_Integration.jl +++ b/test/pde_systems/MOL_1D_Integration.jl @@ -159,3 +159,278 @@ end # @test cumuSumsol ≈ exact atol = 0.36 end + +@testset "Test 03: Partial integro-differential equation with time derivative, with sys transformation" begin + # Equation: ∂u/∂t + u(t,x) + ∫₀ˣ u(t,ξ) dξ = x + x²/2 * (1 - exp(-t - τ)) with τ = 0.2 (alternatively τ = 0) + # On domain: x ∈ [0,2] and t ∈ [0,1] + # Initial condition: u(0,x) = x * (1 - exp(-τ)) + # Boundary condition: u(t,0) = 0 + # Analytical solution: u(t,x) = x * (1 - exp(-t - τ)) + + @parameters t, x + @variables u(..) + Dt = Differential(t) + Dx = Differential(x) + xmin = 0.0 + xmax = 2.0 + + Ix_0_to_x = Integral(x in DomainSets.ClosedInterval(xmin, x)) # basically cumulative sum from 0 to x + + eqs = [Dt(u(t,x)) + u(t,x) + Ix_0_to_x(u(t,x)) ~ x + x^2/2*(1-exp(-t - 0.2))] + bcs = [ + u(0,x) ~ x * (1 - exp(-0.2)), # Initial condition + u(t,0) ~ 0.0] # Boundary condition + domains = [ + t ∈ Interval(0.0, 1.0), + x ∈ Interval(xmin, xmax)] + @named pde_system = PDESystem(eqs, bcs, domains, [t, x], [u(t, x)]) + + # Discretize and compute numerical solution + disc = MOLFiniteDifference([x => 50], t) + prob = discretize(pde_system, disc) + sol = solve(prob, Tsit5(), saveat=0.1) + + # Extract numerical solution at grid points + xdisc = sol[x] + tdisc = sol[t] + usol = sol[u(t,x)] + + # Calculate analytical solution at grid points + analytical_solution(t, x) = x * (1 - exp(-t - 0.2) ) + exact = [analytical_solution(t_, x_) for t_ in tdisc, x_ in xdisc] + + # Compare numerical and analytical solutions + @test usol ≈ exact atol=0.1 + + # # Create animated plot + # using Plots # Add this import + # anim = @animate for (i, t_val) in enumerate(tdisc) + # plot(xdisc, usol[i, :], + # label="Numerical", + # linewidth=2, + # xlabel="x", + # ylabel="u(t,x)", + # title="Solution at t = $(round(t_val, digits=3))", + # ylim=(0, maximum(usol)*1.1), + # xlim=(xmin, xmax)) + + # plot!(xdisc, exact[i, :], + # label="Analytical", + # linewidth=2, + # linestyle=:dash) + + # end + # gif(anim, "integro_diff_solution.gif", fps=10) # Save the animation +end + +@testset "Test 03b: Mirrored version of: Partial integro-differential equation with time derivative, with sys transformation" begin + # Original equation: ∂u/∂t + u(t,x) + ∫₀ˣ u(t,ξ) dξ = x + x²/2 * (1 - exp(-t - τ)) with τ = 0 + # Mirrored version with transformation y = -x, x = -y: + # Equation: ∂u/∂t + u(t,y) + ∫_y ^⁰ u(t,ξ) dξ = -y + y²/2 * (1 - exp(-t)) + # On domain: y ∈ [-2, 0] and t ∈ [0,1] + # Initial condition: u(0,y) = 0 + # Boundary condition: u(t,0) = 0 + # Analytical solution: u(t,y) = -y * (1 - exp(-t)) + + @parameters t, y + @variables u(..) + Dt = Differential(t) + Dy = Differential(y) + ymin = -2.0 + ymax = 0.0 + + # Integral from y to 0 (note the reversed limits) + Iy_to_0 = Integral(y in DomainSets.ClosedInterval(y, ymax)) + + eqs = [Dt(u(t,y)) + u(t,y) + Iy_to_0(u(t,y)) ~ -y + y^2/2*(1-exp(-t))] + bcs = [ + u(0,y) ~ 0.0, # Initial condition + u(t,0) ~ 0.0] # Boundary condition (at y = 0, which corresponds to x = 0) + domains = [ + t ∈ Interval(0.0, 1.0), + y ∈ Interval(ymin, ymax)] + @named pde_system = PDESystem(eqs, bcs, domains, [t, y], [u(t, y)]) + + # Discretize and compute numerical solution + disc = MOLFiniteDifference([y => 50], t) + prob = discretize(pde_system, disc) + sol = solve(prob, Tsit5(), saveat=0.1) + + # Extract numerical solution at grid points + ydisc = sol[y] + tdisc = sol[t] + usol = sol[u(t,y)] + + # Calculate analytical solution at grid points + # Since u(t,x) = x * (1 - exp(-t)) and y = -x, then u(t,y) = -y * (1 - exp(-t)) + analytical_solution(t, y) = -y * (1 - exp(-t)) + exact = [analytical_solution(t_, y_) for t_ in tdisc, y_ in ydisc] + + # Compare numerical and analytical solutions + @test usol ≈ exact atol=0.1 + + # # Create animated plot + # using Plots # Add this import + # anim = @animate for (i, t_val) in enumerate(tdisc) + # plot(ydisc, usol[i, :], + # label="Numerical", + # linewidth=2, + # xlabel="y", + # ylabel="u(t,y)", + # title="Solution at t = $(round(t_val, digits=3))", + # ylim=(0, maximum(usol)*1.1), + # xlim=(ymin, ymax)) + + # plot!(ydisc, exact[i, :], + # label="Analytical", + # linewidth=2, + # linestyle=:dash) + # end + # gif(anim, "integro_diff_solution.gif", fps=10) # Save the animation +end + +@testset "Test 04: PIDE with spatial derivative and source term - domain [0, 4π]" begin + # Equation: ∂u/∂t = ∂u/∂x - u(t,x) - ∫₀ˣ u(t,ξ) dξ + h(x,t) + # where h(x,t) = -sin(x)sin(t) + sin(x)cos(t) - cos(x)cos(t) + cos(t)(1 - cos(x)) + # On domain: x ∈ [0,L]=[0,2π] and t ∈ [0,T]=[0,1] + # Initial condition: u(0,x) = sin(x) + # Boundary conditions: u(t,0) = 0, u(t,2π) = sin(2π)cos(t) = 0 + # Analytical solution: u(t,x) = sin(x)cos(t) + + @parameters t, x + @variables u(..) + Dt = Differential(t) + Dx = Differential(x) + xmin = 0.0 + xmax = 4π + + I0_to_x = Integral(x in DomainSets.ClosedInterval(xmin, x)) # integral from 0 to x + + # Source term h(x,t) that admits the analytical solution + h(x, t) = -sin(x)*sin(t) + sin(x)*cos(t) - cos(x)*cos(t) + cos(t)*(1 - cos(x)) + + # PIDE: ∂u/∂t = ∂u/∂x - u(t,x) - ∫₀ˣ u(t,ξ) dξ + h(x,t) + eqs = [Dt(u(t,x)) ~ Dx(u(t,x)) - u(t,x) - I0_to_x(u(t,x)) + h(x, t)] + bcs = [ + u(0, x) ~ sin(x), # Initial condition + u(t, 0) ~ 0.0, # Boundary condition at x = 0 + u(t, xmax) ~ sin(xmax)*cos(t) # Boundary condition at x = 4π (sin(4π) = 0) + ] + domains = [ + t ∈ Interval(0.0, 1.0), + x ∈ Interval(xmin, xmax)] + @named pde_system = PDESystem(eqs, bcs, domains, [t, x], [u(t, x)]) + + # Discretize and compute numerical solution + disc = MOLFiniteDifference([x => 200], t) + prob = discretize(pde_system, disc) + sol = solve(prob, Tsit5(), saveat=0.1) + + # Extract numerical solution at grid points + xdisc = sol[x] + tdisc = sol[t] + usol = sol[u(t,x)] + + # Calculate analytical solution at grid points + analytical_solution(t, x) = sin(x) * cos(t) + exact = [analytical_solution(t_, x_) for t_ in tdisc, x_ in xdisc] + + # Compare numerical and analytical solutions + @test usol ≈ exact atol=0.36 + + # # Create animated plot + # using Plots # Add this import + # anim = @animate for (i, t_val) in enumerate(tdisc) + # plot(xdisc, usol[i, :], + # label="Numerical", + # linewidth=2, + # xlabel="x", + # ylabel="u(t,x)", + # title="PIDE Solution [0,4π] at t = $(round(t_val, digits=3))", + # ylim=(-1.1, 1.1), + # xlim=(xmin, xmax)) + + # plot!(xdisc, exact[i, :], + # label="Analytical", + # linewidth=2, + # linestyle=:dash) + + # end + # gif(anim, "pide_solution_4pi.gif", fps=10) # Save the animation +end + + +@testset "Test 04b: PIDE with spatial derivative and source term - mirrored domain [-4π, 0]" begin + # Mirrored version with transformation y = -x, so x = -y + # Original equation: ∂u/∂t = ∂u/∂x - u(t,x) - ∫₀ˣ u(t,ξ) dξ + h(x,t) + # Transformed equation: ∂u/∂t = -∂u/∂y - u(t,y) - ∫ₓ⁰ u(t,η) dη + h(-y,t) + # where h(-y,t) = -sin(-y)sin(t) + sin(-y)cos(t) - cos(-y)cos(t) + cos(t)(1 - cos(-y)) + # = sin(y)sin(t) - sin(y)cos(t) - cos(y)cos(t) + cos(t)(1 - cos(y)) + # On domain: y ∈ [-4π, 0] and t ∈ [0,1] + # Initial condition: u(0,y) = sin(-y) = -sin(y) + # Boundary conditions: u(t,0) = 0, u(t,-4π) = sin(-(-4π))cos(t) = sin(4π)cos(t) = 0 + # Analytical solution: u(t,y) = sin(-y)cos(t) = -sin(y)cos(t) + + @parameters t, y + @variables u(..) + Dt = Differential(t) + Dy = Differential(y) + ymin = -4π + ymax = 0.0 + + # Integral from y to 0 (note the reversed limits) + Iy_to_0 = Integral(y in DomainSets.ClosedInterval(y, ymax)) + + # Source term h(-y,t) transformed for y coordinate + h_transformed(y, t) = sin(y)*sin(t) - sin(y)*cos(t) - cos(y)*cos(t) + cos(t)*(1 - cos(y)) + + # PIDE: ∂u/∂t = -∂u/∂y - u(t,y) - ∫ʸ⁰ u(t,η) dη + h(-y,t) + eqs = [Dt(u(t,y)) ~ -Dy(u(t,y)) - u(t,y) - Iy_to_0(u(t,y)) + h_transformed(y, t)] + bcs = [ + u(0, y) ~ -sin(y), # Initial condition: u(0,y) = sin(-y) = -sin(y) + u(t, 0) ~ 0.0, # Boundary condition at y = 0 + u(t, ymin) ~ -sin(ymin)*cos(t) # Boundary condition at y = -4π (sin(-4π) = 0) + ] + domains = [ + t ∈ Interval(0.0, 1.0), + y ∈ Interval(ymin, ymax)] + @named pde_system = PDESystem(eqs, bcs, domains, [t, y], [u(t, y)]) + + # Discretize and compute numerical solution + disc = MOLFiniteDifference([y => 200], t) + prob = discretize(pde_system, disc) + sol = solve(prob, Tsit5(), saveat=0.1) + + # Extract numerical solution at grid points + ydisc = sol[y] + tdisc = sol[t] + usol = sol[u(t,y)] + + # Calculate analytical solution at grid points + # Analytical solution: u(t,y) = sin(-y)cos(t) = -sin(y)cos(t) + analytical_solution(t, y) = -sin(y) * cos(t) + exact = [analytical_solution(t_, y_) for t_ in tdisc, y_ in ydisc] + + # Compare numerical and analytical solutions + @test usol ≈ exact atol=0.36 + + # # Create animated plot + # using Plots # Add this import + # anim = @animate for (i, t_val) in enumerate(tdisc) + # plot(ydisc, usol[i, :], + # label="Numerical", + # linewidth=2, + # xlabel="y", + # ylabel="u(t,y)", + # title="PIDE Solution [-4π,0] at t = $(round(t_val, digits=3))", + # ylim=(-1.1, 1.1), + # xlim=(ymin, ymax)) + + # plot!(ydisc, exact[i, :], + # label="Analytical", + # linewidth=2, + # linestyle=:dash) + + # end + # gif(anim, "pide_solution_mirrored_4pi.gif", fps=10) # Save the animation +end