Merge pull request #483 from SciML/fixmaster

fix master

Author: ChrisRackauckas

src/MOLFiniteDifference/MOL_discretization.jl

@@ -334,9 +334,24 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
                                             Num(substitute(substitute(*(~~a..., ~~b...), r_mid_dep(II, j, k, 1)), r_mid_indep(II, j, 1)))],
                                             [-b1(II, j, k), b2(II, j, k)])
                                         for (j, iv) in enumerate(nottime) for (k, dv) in enumerate(depvars)]
-                spherical_deriv_rules = [@rule (~a) * (iv^-2) * ($(Differential(iv))((iv^2)*$(Differential(iv))(dv))) =>
-                                         (~a) * central_deriv_spherical(II, j, k)
+
+                cartesian_deriv_rules = vcat(vec(cartesian_deriv_rules),vec(
+                                        [@rule ($(Differential(iv))($(Differential(iv))(dv)/~a)) =>
+                                        dot([Num(substitute(substitute(1/~a, r_mid_dep(II, j, k, -1)), r_mid_indep(II, j, -1))),
+                                            Num(substitute(substitute(1/~a, r_mid_dep(II, j, k, 1)), r_mid_indep(II, j, 1)))],
+                                            [-b1(II, j, k), b2(II, j, k)])
+                                        for (j, iv) in enumerate(nottime) for (k, dv) in enumerate(depvars)]))
+
+                spherical_deriv_rules = [@rule *(~~a, ($(Differential(iv))((iv^2)*$(Differential(iv))(dv))), ~~b) / (iv^2) =>
+                                         *(~a..., central_deriv_spherical(II, j, k), ~b...)
                                          for (j, iv) in enumerate(nottime) for (k, dv) in enumerate(depvars)]
+
+                # r^-2 needs to be handled separately
+                spherical_deriv_rules = vcat(vec(spherical_deriv_rules),vec(
+                                        [@rule *(~~a, (iv^-2) * ($(Differential(iv))((iv^2)*$(Differential(iv))(dv))), ~~b) =>
+                                         *(~a..., central_deriv_spherical(II, j, k), ~b...)
+                                         for (j, iv) in enumerate(nottime) for (k, dv) in enumerate(depvars)]))
+
                 rhs_arg = istree(eq.rhs) && (SymbolicUtils.operation(eq.rhs) == +) ? SymbolicUtils.arguments(eq.rhs) : [eq.rhs]
                 lhs_arg = istree(eq.lhs) && (SymbolicUtils.operation(eq.lhs) == +) ? SymbolicUtils.arguments(eq.lhs) : [eq.lhs]
                 nonlinlap_rules = []