From 13b0e6baecf6fe824f1ff60ed4502cda613a4308 Mon Sep 17 00:00:00 2001 From: Ian Hinder Date: Sat, 17 Sep 2011 21:41:40 +0200 Subject: McLachlan_BSSN.m: Make finite difference order a run-time parameter --- m/McLachlan_BSSN.m | 188 +++++++++++++---------------------------------------- 1 file changed, 46 insertions(+), 142 deletions(-) (limited to 'm') diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m index a999fb6..cc5ef93 100644 --- a/m/McLachlan_BSSN.m +++ b/m/McLachlan_BSSN.m @@ -33,156 +33,55 @@ KD = KroneckerDelta; derivatives = { - PDstandardNth[i_] -> StandardCenteredDifferenceOperator[1,derivOrder/2,i], - PDstandardNth[i_,i_] -> StandardCenteredDifferenceOperator[2,derivOrder/2,i], - PDstandardNth[i_,j_] -> StandardCenteredDifferenceOperator[1,derivOrder/2,i] * - StandardCenteredDifferenceOperator[1,derivOrder/2,j], + PDstandardNth[i_] -> StandardCenteredDifferenceOperator[1,fdOrder/2,i], + PDstandardNth[i_,i_] -> StandardCenteredDifferenceOperator[2,fdOrder/2,i], + PDstandardNth[i_,j_] -> StandardCenteredDifferenceOperator[1,fdOrder/2,i] * + StandardCenteredDifferenceOperator[1,fdOrder/2,j], PDdissipationNth[i_] -> - spacing[i]^(derivOrder+1) / 2^(derivOrder+2) * - StandardCenteredDifferenceOperator[derivOrder+2,derivOrder/2+1,i], + spacing[i]^(fdOrder+1) / 2^(fdOrder+2) * + StandardCenteredDifferenceOperator[fdOrder+2,fdOrder/2+1,i], (* PD: These come from my mathematica notebook "Upwind-Kranc-Convert.nb" that converts upwinding finite differencing operators generated by StandardUpwindDifferenceOperator into this form *) - Switch[derivOrder, - 2, - PDupwindNth[1] -> (dir[1]*(-3 + 4*shift[1]^dir[1] - - shift[1]^(2*dir[1])))/(2*spacing[1]), - 4, - PDupwindNth[1] -> (dir[1]*(-10 - 3/shift[1]^dir[1] + 18*shift[1]^dir[1] - - 6*shift[1]^(2*dir[1]) + - shift[1]^(3*dir[1])))/(12*spacing[1]), - 6, - PDupwindNth[1] -> (dir[1]*(-35 + 2/shift[1]^(2*dir[1]) - - 24/shift[1]^dir[1] + 80*shift[1]^dir[1] - - 30*shift[1]^(2*dir[1]) + 8*shift[1]^(3*dir[1]) - - shift[1]^(4*dir[1])))/(60*spacing[1]), - 8, - PDupwindNth[1] -> (dir[1]*(-378 - 5/shift[1]^(3*dir[1]) + - 60/shift[1]^(2*dir[1]) - 420/shift[1]^dir[1] + - 1050*shift[1]^dir[1] - 420*shift[1]^(2*dir[1]) + - 140*shift[1]^(3*dir[1]) - 30*shift[1]^(4*dir[1]) + - 3*shift[1]^(5*dir[1])))/(840*spacing[1])], - Switch[derivOrder, - 2, - PDupwindNth[2] -> (dir[2]*(-3 + 4*shift[2]^dir[2] - - shift[2]^(2*dir[2])))/(2*spacing[2]), - 4, - PDupwindNth[2] -> (dir[2]*(-10 - 3/shift[2]^dir[2] + 18*shift[2]^dir[2] - - 6*shift[2]^(2*dir[2]) + - shift[2]^(3*dir[2])))/(12*spacing[2]), - 6, - PDupwindNth[2] -> (dir[2]*(-35 + 2/shift[2]^(2*dir[2]) - - 24/shift[2]^dir[2] + 80*shift[2]^dir[2] - - 30*shift[2]^(2*dir[2]) + 8*shift[2]^(3*dir[2]) - - shift[2]^(4*dir[2])))/(60*spacing[2]), - 8, - PDupwindNth[2] -> (dir[2]*(-378 - 5/shift[2]^(3*dir[2]) + - 60/shift[2]^(2*dir[2]) - 420/shift[2]^dir[2] + - 1050*shift[2]^dir[2] - 420*shift[2]^(2*dir[2]) + - 140*shift[2]^(3*dir[2]) - 30*shift[2]^(4*dir[2]) + - 3*shift[2]^(5*dir[2])))/(840*spacing[2])], - Switch[derivOrder, - 2, - PDupwindNth[3] -> (dir[3]*(-3 + 4*shift[3]^dir[3] - - shift[3]^(2*dir[3])))/(2*spacing[3]), - 4, - PDupwindNth[3] -> (dir[3]*(-10 - 3/shift[3]^dir[3] + 18*shift[3]^dir[3] - - 6*shift[3]^(2*dir[3]) + - shift[3]^(3*dir[3])))/(12*spacing[3]), - 6, - PDupwindNth[3] -> (dir[3]*(-35 + 2/shift[3]^(2*dir[3]) - - 24/shift[3]^dir[3] + 80*shift[3]^dir[3] - - 30*shift[3]^(2*dir[3]) + 8*shift[3]^(3*dir[3]) - - shift[3]^(4*dir[3])))/(60*spacing[3]), - 8, - PDupwindNth[3] -> (dir[3]*(-378 - 5/shift[3]^(3*dir[3]) + - 60/shift[3]^(2*dir[3]) - 420/shift[3]^dir[3] + - 1050*shift[3]^dir[3] - 420*shift[3]^(2*dir[3]) + - 140*shift[3]^(3*dir[3]) - 30*shift[3]^(4*dir[3]) + - 3*shift[3]^(5*dir[3])))/(840*spacing[3])], - - (* TODO: make these higher order stencils *) - PDonesided[1] -> dir[1] (-1 + shift[1]^dir[1]) / spacing[1], - PDonesided[2] -> dir[2] (-1 + shift[2]^dir[2]) / spacing[2], - PDonesided[3] -> dir[3] (-1 + shift[3]^dir[3]) / spacing[3] + Sequence@@Flatten[Table[ + {PDupwindNth[i] -> Switch[fdOrder, + 2, (dir[i]*(-3 + 4*shift[i]^dir[i] - shift[i]^(2*dir[i])))/(2*spacing[i]), + 4, (dir[i]*(-10 - 3/shift[i]^dir[i] + 18*shift[i]^dir[i] - + 6*shift[i]^(2*dir[i]) + shift[i]^(3*dir[i])))/(12*spacing[i]), + 6, (dir[i]*(-35 + 2/shift[i]^(2*dir[i]) - 24/shift[i]^dir[i] + 80*shift[i]^dir[i] - + 30*shift[i]^(2*dir[i]) + 8*shift[i]^(3*dir[i]) - shift[i]^(4*dir[i])))/(60*spacing[i]), + 8, (dir[i]*(-378 - 5/shift[i]^(3*dir[i]) + 60/shift[i]^(2*dir[i]) - 420/shift[i]^dir[i] + + 1050*shift[i]^dir[i] - 420*shift[i]^(2*dir[i]) + 140*shift[i]^(3*dir[i]) - 30*shift[i]^(4*dir[i]) + + 3*shift[i]^(5*dir[i])))/(840*spacing[i])], + + PDupwindNthAnti[i] -> Switch[fdOrder, + 2, (+1 shift[i]^(-2) -4 shift[i]^(-1) +0 shift[i]^( 0) +4 shift[i]^(+1) -1 shift[i]^(+2)) / (4 spacing[i]), + 4, (-1 shift[i]^(-3) +6 shift[i]^(-2) -21 shift[i]^(-1 )+0 shift[i]^( 0) +21 shift[i]^(+1) + -6 shift[i]^(+2) +1 shift[i]^(+3)) / (24 spacing[i]), + 6, (+1 shift[i]^(-4) -8 shift[i]^(-3) +32 shift[i]^(-2) -104 shift[i]^(-1) +0 shift[i]^( 0) + +104 shift[i]^(+1) -32 shift[i]^(+2) +8 shift[i]^(+3) -1 shift[i]^(+4)) / (120 spacing[i]), + 8, (-3 shift[i]^(-5) +30 shift[i]^(-4) -145 shift[i]^(-3) +480 shift[i]^(-2) -1470 shift[i]^(-1) + +0 shift[i]^( 0) +1470 shift[i]^(+1) -480 shift[i]^(+2) +145 shift[i]^(+3) -30 shift[i]^(+4) + +3 shift[i]^(+5)) / (1680 spacing[i])], + + PDupwindNthSymm[i] -> Switch[fdOrder, + 2, (-1 shift[i]^(-2) +4 shift[i]^(-1) -6 shift[i]^( 0) +4 shift[i]^(+1) -1 shift[i]^(+2)) / (4 spacing[i]), + 4, (+1 shift[i]^(-3) -6 shift[i]^(-2) +15 shift[i]^(-1) -20 shift[i]^( 0) +15 shift[i]^(+1) + -6 shift[i]^(+2) +1 shift[i]^(+3)) / (24 spacing[i]), + 6, (-1 shift[i]^(-4) +8 shift[i]^(-3) - 28 shift[i]^(-2)+56 shift[i]^(-1)-70 shift[i]^( 0) + +56 shift[i]^(+1) -28 shift[i]^(+2) +8 shift[i]^(+3) -1 shift[i]^(+4)) / (120 spacing[i]), + 8, (+3 shift[i]^(-5) -30 shift[i]^(-4) +135 shift[i]^(-3) -360 shift[i]^(-2) +630 shift[i]^(-1) + -756 shift[i]^( 0) +630 shift[i]^(+1) -360 shift[i]^(+2) +135 shift[i]^(+3) -30 shift[i]^(+4) + +3 shift[i]^(+5)) / (1680 spacing[i])], + + (* TODO: make these higher order stencils *) + PDonesided[i] -> dir[i] (-1 + shift[i]^dir[i]) / spacing[i]} /. i->j, {j,1,3}],1] }; -derivatives = Join[derivatives, - Flatten[Table[ - Switch[derivOrder, - 2, - {PDupwindNthAnti[i] -> ( +1 shift[i]^(-2) - -4 shift[i]^(-1) - +0 shift[i]^( 0) - +4 shift[i]^(+1) - -1 shift[i]^(+2)) / (4 spacing[i]), - PDupwindNthSymm[i] -> ( -1 shift[i]^(-2) - +4 shift[i]^(-1) - -6 shift[i]^( 0) - +4 shift[i]^(+1) - -1 shift[i]^(+2)) / (4 spacing[i])}, - 4, - {PDupwindNthAnti[i] -> ( -1 shift[i]^(-3) - +6 shift[i]^(-2) - -21 shift[i]^(-1) - +0 shift[i]^( 0) - +21 shift[i]^(+1) - -6 shift[i]^(+2) - +1 shift[i]^(+3)) / (24 spacing[i]), - PDupwindNthSymm[i] -> ( +1 shift[i]^(-3) - -6 shift[i]^(-2) - +15 shift[i]^(-1) - -20 shift[i]^( 0) - +15 shift[i]^(+1) - -6 shift[i]^(+2) - +1 shift[i]^(+3)) / (24 spacing[i])}, - 6, - {PDupwindNthAnti[i] -> ( +1 shift[i]^(-4) - -8 shift[i]^(-3) - +32 shift[i]^(-2) - -104 shift[i]^(-1) - +0 shift[i]^( 0) - +104 shift[i]^(+1) - -32 shift[i]^(+2) - +8 shift[i]^(+3) - -1 shift[i]^(+4)) / (120 spacing[i]), - PDupwindNthSymm[i] -> ( -1 shift[i]^(-4) - +8 shift[i]^(-3) - -28 shift[i]^(-2) - +56 shift[i]^(-1) - -70 shift[i]^( 0) - +56 shift[i]^(+1) - -28 shift[i]^(+2) - +8 shift[i]^(+3) - -1 shift[i]^(+4)) / (120 spacing[i])}, - 8, - {PDupwindNthAnti[i] -> ( -3 shift[i]^(-5) - +30 shift[i]^(-4) - -145 shift[i]^(-3) - +480 shift[i]^(-2) - -1470 shift[i]^(-1) - +0 shift[i]^( 0) - +1470 shift[i]^(+1) - -480 shift[i]^(+2) - +145 shift[i]^(+3) - -30 shift[i]^(+4) - +3 shift[i]^(+5)) / (1680 spacing[i]), - PDupwindNthSymm[i] -> ( +3 shift[i]^(-5) - -30 shift[i]^(-4) - +135 shift[i]^(-3) - -360 shift[i]^(-2) - +630 shift[i]^(-1) - -756 shift[i]^( 0) - +630 shift[i]^(+1) - -360 shift[i]^(+2) - +135 shift[i]^(+3) - -30 shift[i]^(+4) - +3 shift[i]^(+5)) / (1680 spacing[i])}], - {i,1,3}]] -]; +(* Put[derivatives, "derivatives.m"]; *) PD = PDstandardNth; PDu = PDupwindNth; @@ -1242,6 +1141,11 @@ intParameters = AllowedValues -> {{Value -> "0", Description -> "phi method"}, {Value -> "1", Description -> "W method"}}, Default -> 0 + }, + { + Name -> fdOrder, + Default -> derivOrder, + AllowedValues -> {2,4,6,8} } }; @@ -1346,7 +1250,7 @@ CreateKrancThornTT [groups, ".", BSSN, PartialDerivatives -> derivatives, EvolutionTimelevels -> evolutionTimelevels, DefaultEvolutionTimelevels -> 3, - UseJacobian -> useJacobian, + UseJacobian -> True, UseLoopControl -> True, UseVectors -> True, InheritedImplementations -> inheritedImplementations, -- cgit v1.2.3