stash

Author: xtalax

src/MOLFiniteDifference.jl

@@ -9,7 +9,7 @@ end
 function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, grid_align=CenterAlignedGrid())
 
     if approx_order % 2 != 0
-        warn("Discretization approx_order must be even, rounding up to $(approx_order+1)")
+        @warn "Discretization approx_order must be even, rounding up to $(approx_order+1)"
     end
     return MOLFiniteDifference{typeof(grid_align)}(dxs, time, approx_order, grid_align)
 end

src/bcs/generate_bc_eqs.jl

@@ -94,7 +94,6 @@ function BoundaryHandler!!(u0, bceqs, bcs, s::DiscreteSpace, depvar_ops, tspan,
                 end
 
                 @assert boundary !== nothing "Boundary condition $bc is not on a boundary of the domain, or is not a valid boundary condition"
-                @show bc
                 push!(bceqs, vec(map(s.Iedge[x_][idx(boundary)]) do II
                     rules = generate_bc_rules(II, derivweights, s, bc, u_, x_, boundary)
 
@@ -123,7 +122,6 @@ function generate_bc_rules(II, derivweights, s::DiscreteSpace{N,M,G}, bc, u_, x_
             break
         end
     end
-    @show depvarderivbcmaps
     fd_rules = generate_finite_difference_rules(II, s, bc, derivweights)
     varrules = axiesvals(s, x_, II)
 

src/bcs/scratch.jl

@@ -4,57 +4,65 @@ using DiffEqOperators
 """
 Inner function for `get_half_offset_weights_and_stencil` (see below).
 """
-function _get_weights_and_stencil(D, I)
+function _get_weights_and_stencil(D, I, len=0)
+    if len == 20
+        clip = false
+    else
+        clip = true
+    end
     # k => i of itap - 
     # offset is important due to boundary proximity
     # The low boundary coeffs has a heirarchy of coefficients following: number of indices from boundary -> which half offset point does it correspond to -> weights
     if I <= (D.boundary_point_count)
         println("low boundary")
         weights = D.low_boundary_coefs[I]    
+        weight = I
         offset = 1 - I
         Itap = [I + (i+offset) for i in 0:(D.boundary_stencil_length-1)]
-    elseif I > (20 - D.boundary_point_count)
+    elseif I > (len - D.boundary_point_count)
         println("high boundary")
 
-        weights = D.high_boundary_coefs[20-I+1]
-        offset = 20 - I
-        Itap = [I + (i+offset) for i in (-D.boundary_stencil_length+1):1:0]
+        weights = D.high_boundary_coefs[len-I + !clip]
+        weight = -(len - I + !clip)
+        offset = len - I
+        Itap = [I + (i+offset+clip) for i in (-D.boundary_stencil_length+1):1:0]
     else
         println("interior")
-
+        weight = 0
         weights = D.stencil_coefs
         Itap = [I + i for i in (1-div(D.stencil_length,2)):(div(D.stencil_length,2))]
     end    
 
-    return (weights, Itap)
+    return (weights, Itap, weight)
 end
 
 
-function clip(I, n)
-    return I>(n/2) ? I-1 : I+1
+function clip(I)
+    return I-1 
 
 end
 
 for I in [2,19, 3, 18, 10]
 
-    D =  CompleteHalfCenteredDifference(1, 6, 1.0)
+    D =  CompleteHalfCenteredDifference(1, 2, 1.0)
 
-     
+    x = 1:20
 
-    outerweights, outerstencil = _get_weights_and_stencil(D, clip(I,20))
+    outerweights, outerstencil, outerweight = _get_weights_and_stencil(D, clip(I), 19)
 
     # Get the correct weights and stencils for this II
 
     @show I
+    @show outerweight
     @show outerweights
     @show outerstencil
 
-    
     for II in outerstencil
         println()
         println("II: ", II)
         println()
-        weights, stencil = _get_weights_and_stencil(D, II)
+        weights, stencil, weight = _get_weights_and_stencil(D, II, 20)
+        @show weight
         println(weights)
         println(stencil)
     end

src/bcs/scratch2.jl

@@ -0,0 +1,21 @@
+using ModelingToolkit, MethodOfLines
+
+@parameters t x
+@variables u(..)
+Dx = Differential(x)
+Dt = Differential(t)
+t_min= 0.
+t_max = 2.
+x_min = 0.
+x_max = 2.
+c = 1.0
+a = 1.0
+
+# Equation
+eq = Dt(u(t,x)) ~ Dx(u(t,x) * Dx(u(t,x)))
+
+terms = MethodOfLines.split_additive_terms(eq)
+
+r = @rule ($(Differential(x))(u(t,x) * $(Differential(x))(u(t,x)))) => 0
+
+@show r(terms[2])

src/discretization/differential_discretizer.jl

@@ -71,8 +71,15 @@ Get the weights and stencil for the inner half offset centered difference for th
 Does not discretize so that the weights can be used in a replacement rule
 TODO: consider refactoring this to harmonize with centered difference
 """
-function get_half_offset_weights_and_stencil(D, II, s, jx)
+function get_half_offset_weights_and_stencil(D::DerivativeOperator, II::CartesianIndex, s::DiscreteSpace, jx, len::Int = 0)
     j, x = jx
+    # Check if the index is clipped because we need the outerweights
+    if len == 0
+        len = length(s, x)
+        clip = false
+    else
+        clip = true
+    end
     # unit index in direction of the derivative
     I1 = unitindices(nparams(s))[j] 
     # offset is important due to boundary proximity
@@ -81,10 +88,12 @@ function get_half_offset_weights_and_stencil(D, II, s, jx)
         weights = D.low_boundary_coefs[II[j]]    
         offset = 1 - II[j]
         Itap = [II + (i+offset)*I1 for i in 0:(D.boundary_stencil_length-1)]
-    elseif II[j] > (length(s, x) - D.boundary_point_count)
-        weights = D.high_boundary_coefs[length(s, x)-II[j]+1]
-        offset = length(s, x) - II[j]
-        Itap = [II + (i+offset)*I1 for i in (-D.boundary_stencil_length+1):1:0]
+    elseif II[j] > (len - D.boundary_point_count)
+        # we need !clip to shift the index, as unclipped indices start from 1 not 2
+        weights = D.high_boundary_coefs[len-II[j] + !clip]
+        offset = len - II[j]
+        #use clip to shift the stencil to the right place
+        Itap = [II + (i+offset+clip)*I1 for i in (-D.boundary_stencil_length+1):1:0]
     else
         weights = D.stencil_coefs
         Itap = [II + i*I1 for i in (1-div(D.stencil_length,2)):(div(D.stencil_length,2))]

src/discretization/generate_finite_difference_rules.jl

@@ -4,7 +4,12 @@
 Interpolate gridpoints by taking the average of the values of the discrete points, or if the offset is outside the grid, extrapolate the value with dx.
 """
 @inline function interpolate_discrete_param(i, s, itap, x, bpc)
-    return s.grid[x][i+itap]+s.dxs[x]*.5
+    if i+itap > (length(s, x) - bpc)
+        return s.grid[x][i+itap]-s.dxs[x]*.5        
+
+    else
+        return s.grid[x][i+itap]+s.dxs[x]*.5        
+    end   
 end
 
 """
@@ -46,7 +51,7 @@ function cartesian_nonlinear_laplacian(expr, II, derivweights, s, x, u)
     inner_interpolater = derivweights.interpmap[x]
 
     # Get the outer weights and stencil. clip() essentially removes a point from either end of the grid, for this reason this function is only defined on the interior, not in bcs
-    outerweights, outerstencil = get_half_offset_weights_and_stencil(D_inner, clip(II, s, j, D_inner.boundary_point_count), s, jx)
+    outerweights, outerstencil = get_half_offset_weights_and_stencil(D_inner, clip(II, s, j, D_inner.boundary_point_count), s, jx, length(s, x) - 1)
 
     # Get the correct weights and stencils for this II
     inner_deriv_weights_and_stencil = [get_half_offset_weights_and_stencil(D_inner, I, s, jx) for I in outerstencil]

test/components/finite_diff_schemes.jl

@@ -80,7 +80,7 @@ domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)]
 
     s = MethodOfLines.DiscreteSpace(domains, depvars, indvars, x̄, disc)
 
-    derivweights = MethodOfLines.DifferentialDiscretizer(pde, s, disc)
+    derivweights = MethodOfLines.DifferentialDiscretizer(pde, bcs, s, disc)
 
     II = s.Igrid[10]
 
@@ -118,7 +118,7 @@ end
 
     s = MethodOfLines.DiscreteSpace(domains, depvars, indvars, x̄, disc)
 
-    derivweights = MethodOfLines.DifferentialDiscretizer(pde, s, disc)
+    derivweights = MethodOfLines.DifferentialDiscretizer(pde, bcs, s, disc)
 
     for II in s.Igrid[2:end-1]
         #TODO Test Interpolation of params