From 077267b146ca7e73522c671d3c6892a1bd1cbe7c Mon Sep 17 00:00:00 2001 From: Ian Hinder Date: Wed, 24 Mar 2010 01:14:35 +0100 Subject: Add wave equation example script and parameter file --- Examples/Wave.m | 216 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 216 insertions(+) create mode 100644 Examples/Wave.m (limited to 'Examples/Wave.m') diff --git a/Examples/Wave.m b/Examples/Wave.m new file mode 100644 index 0000000..0c19475 --- /dev/null +++ b/Examples/Wave.m @@ -0,0 +1,216 @@ + +Get["KrancThorn`"]; +SetEnhancedTimes[False]; + +(********************************************************) +(* Derivatives *) +(********************************************************) + +derivatives = +{ + PDstandard2nd[i_] -> StandardCenteredDifferenceOperator[1,1,i], + PDstandard2nd[i_, i_] -> StandardCenteredDifferenceOperator[2,1,i], + PDstandard2nd[i_, j_] -> StandardCenteredDifferenceOperator[1,1,i] StandardCenteredDifferenceOperator[1,1,j], + + PDstandard4th[i_] -> StandardCenteredDifferenceOperator[1,2,i], + PDstandard4th[i_, i_] -> StandardCenteredDifferenceOperator[2,2,i], + PDstandard4th[i_, j_] -> StandardCenteredDifferenceOperator[1,2,i] StandardCenteredDifferenceOperator[1,2,j], + + PDonesided2nd[1] -> dir[1] (-shift[1]^(2 dir[1]) + 4 shift[1]^dir[1] - 3 )/(2 spacing[1]), + PDonesided2nd[2] -> dir[2] (-shift[2]^(2 dir[2]) + 4 shift[2]^dir[2] - 3 )/(2 spacing[2]), + PDonesided2nd[3] -> dir[3] (-shift[3]^(2 dir[3]) + 4 shift[3]^dir[3] - 3 )/(2 spacing[3]), + + Diss2nd[] -> - diss Sum[spacing[i]^3 (DPlus[i] DMinus[i])^2, {i, 1, 3}], + Diss4th[] -> diss Sum[spacing[i]^5 (DPlus[i] DMinus[i])^3, {i, 1, 3}], + + PDzero[i_] -> DZero[i], + PDzero[i_, j_] -> DZero[i] DZero[j], + + PDplus[i_] -> DPlus[i] +}; + +Map[DefineTensor, {norm, dir}]; + +zerodims = {}; + +PD = PDstandard2nd; +PDnorm = PDplus; + +evolvedGroup = {"evolved", {phi, pi}}; +normGroup = {"norms", {VL2, VDP, EL2}}; +exactGroup = {"exact", {phiExact, piExact}}; +errorGroup = {"errors", {phiError, piError}}; + +groups = {evolvedGroup, normGroup, exactGroup, errorGroup}; + +wp = {Name -> periodicity, Default -> 1}; +wa = {Name -> amplitude, Default -> 1}; + +nNorm = Sqrt[n1^2+n2^2+n3^2]; +nX = n1 x + n2 y + n3 z; + +(* Sine wave exact solution *) + +exactSineCalc = +{ + Name -> "wave_exact_sine", + Before -> {"import_exact"}, + ConditionalOnKeyword -> {"initial_data", "sine"}, + Schedule -> {"AT INITIAL before import_exact","AT POSTSTEP before calc_errors"}, + Shorthands -> {piconst}, + Equations -> + { + piconst -> N[Pi,20], + phiExact -> amplitude Sin[2 piconst / periodicity (nX - nNorm t)], + piExact -> -2 piconst / periodicity nNorm amplitude Cos[2 piconst / periodicity (nX - nNorm t)] + } +} + +(* Radial profile of exact Gaussian solution *) +f[x_] := x^3 Exp[-x^2/nSigma^2]; + +exactGaussianCalc = +{ + Name -> "wave_exact_gaussian", + Before -> {"import_exact"}, + ConditionalOnKeyword -> {"initial_data", "gaussian"}, + Schedule -> {"AT INITIAL before import_exact","AT POSTSTEP before calc_errors"}, + Shorthands -> {piconst,rEps}, + Equations -> + { + rEps -> (r^4+(10^(-6))^4)^(1/4), + + phiExact -> (f[t+t0+r] - f[t+t0-r]) / rEps, + piExact -> (D[f[t+t0+r],t] - D[f[t+t0-r],t]) / rEps + } +} + +(* Import the exact solution as the initial data *) + +importerEquations = +{ + phi -> phiExact, + pi -> piExact +}; + +importerCalc = +{ + Name -> "wave_import_exact", + Schedule -> {"at INITIAL as import_exact"}, + Equations -> importerEquations +} + +(* The evolution equations *) + +evolveEquations = +{ + dot[phi] -> pi, + dot[pi] -> PD[phi,li,ui] +} + +evolveCalc = +{ + Name -> "wave_evolve", + Schedule -> {"in MoL_CalcRHS as evolve"}, + Where -> Interior, + Equations -> evolveEquations +} + +(* Evaluate the errors *) + +errorEquations = +{ + phiError -> phi - phiExact, + piError -> pi - piExact +} + +errorCalc = +{ + Name -> "wave_calc_errors", + Schedule -> {"at ANALYSIS as calc_errors"}, + Equations -> errorEquations +} + +(* Evaluate the norms *) + +normEquations = +{ + VL2squared -> phi^2 + pi^2, + VL2 -> Sqrt[VL2squared], + VDPsquared -> pi^2 + PDnorm[phi,li] PDnorm[phi,ui], + VDP -> Sqrt[VDPsquared], + EL2squared -> phiError^2 + piError^2, + EL2 -> Sqrt[EL2squared] +} + +normCalc = +{ + Name -> "wave_calc_norm", + Where -> Interior, + Schedule -> {"at ANALYSIS as calc_norm"}, + Shorthands -> {VL2squared, VDPsquared, EL2squared}, + Equations -> normEquations +} + +(**************************************************************************************) +(* Boundary conditions *) +(**************************************************************************************) + +boundaryParam = +{ + Name -> "boundary_condition", + Default -> "radiative", + AllowedValues -> {"none", "radiative"} +}; + +boundaryCalc = +{ + Name -> "wave_boundary", + Schedule -> {"in MoL_RHSBoundaries"}, + ConditionalOnKeyword -> {"boundary_condition", "radiative"}, + Where -> Boundary, + Shorthands -> {norm[ui], dir[ui]}, + Equations -> + { + norm1 -> -x/r, + norm2 -> -y/r, + norm3 -> -z/r, + + dir[ui] -> Sign[norm[ui]], + + (* \partial_t u = - (u - u_0) / r - \partial_r u + for u_0 some background solution. In this case, Minkowski in Cartesian + coordinates. *) + + dot[phi] -> -(phi - 0) / r + norm[ui] PDonesided2nd[phi, li], + dot[pi] -> -(pi - 0) / r + norm[ui] PDonesided2nd[pi, li] + } +}; + +(* Construct the thorn *) + +calculations = {exactSineCalc, exactGaussianCalc, importerCalc, evolveCalc, errorCalc, normCalc, boundaryCalc}; + +declaredGroups = Map[groupName, {evolvedGroup, exactGroup, errorGroup, + normGroup}]; + +idParam = +{ + Name -> "initial_data", + Default -> "gaussian", + AllowedValues -> {"gaussian", "sine"} +}; + +keywordParameters = +{ + idParam, + boundaryParam +}; + +CreateKrancThornTT[groups, ".", "Wave", + Calculations -> calculations, + DeclaredGroups -> declaredGroups, + PartialDerivatives -> derivatives, + ZeroDimensions -> zerodims, + KeywordParameters -> keywordParameters, + RealParameters -> {wp, wa, {Name -> n1, Default -> 1}, n2, n3, nSigma, r0, t0, x0}] -- cgit v1.2.3