aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Examples/Wave.m216
-rw-r--r--Examples/wave.par67
2 files changed, 283 insertions, 0 deletions
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}]
diff --git a/Examples/wave.par b/Examples/wave.par
new file mode 100644
index 0000000..e28a05c
--- /dev/null
+++ b/Examples/wave.par
@@ -0,0 +1,67 @@
+
+Cactus::cctk_final_time = 1
+Cactus::terminate = "time"
+
+ActiveThorns = "IOUtil Carpet CarpetLib CarpetSlab CoordBase CoordBase SymBase CartGrid3D Slab CarpetIOBasic CarpetIOASCII Wave Time MoL Periodic Boundary GenericFD CarpetReduce LoopControl"
+
+CoordBase::domainsize = minmax
+
+CoordBase::boundary_size_x_lower = 1
+CoordBase::boundary_size_y_lower = 1
+CoordBase::boundary_size_z_lower = 1
+CoordBase::boundary_shiftout_x_lower = 1
+CoordBase::boundary_shiftout_y_lower = 1
+CoordBase::boundary_shiftout_z_lower = 1
+
+CoordBase::boundary_size_x_upper = 1
+CoordBase::boundary_size_y_upper = 1
+CoordBase::boundary_size_z_upper = 1
+CoordBase::boundary_shiftout_x_upper = 0
+CoordBase::boundary_shiftout_y_upper = 0
+CoordBase::boundary_shiftout_z_upper = 0
+
+CartGrid3D::type = "coordbase"
+CartGrid3D::domain = "full"
+CartGrid3D::avoid_origin = "no"
+
+CoordBase::xmin = 0
+CoordBase::ymin = 0
+CoordBase::zmin = 0
+CoordBase::xmax = 1
+CoordBase::ymax = 1
+CoordBase::zmax = 1
+CoordBase::dx = 0.05
+CoordBase::dy = 0.05
+CoordBase::dz = 0.05
+
+driver::ghost_size = 1
+Carpet::domain_from_coordbase = "yes"
+Carpet::poison_new_timelevels = "yes"
+Carpet::check_for_poison = "no"
+Carpet::poison_value = 113
+
+Carpet::max_refinement_levels = 1
+Carpet::prolongation_order_space = 3
+Carpet::prolongation_order_time = 2
+Carpet::init_each_timelevel = no
+Carpet::print_timestats_every = 0
+
+Periodic::periodic = yes
+
+Time::dtfac = 0.5
+
+MethodOfLines::ode_method = "RK4"
+MethodOfLines::MoL_NaN_Check = "no"
+MethodOfLines::initial_data_is_crap = "yes"
+MethodOfLines::MoL_Intermediate_Steps = 4
+MethodOfLines::MoL_Num_Scratch_Levels = 1
+
+IO::out_dir = $parfile
+
+CarpetIOBasic::outInfo_every = 1
+CarpetIOBasic::outInfo_vars = ""
+
+IOASCII::out1D_every = 1
+IOASCII::out1D_vars = "Wave::phi Wave::phiExact"
+Wave::initial_data = "sine"
+Wave::evolved_bound = "none"
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-16 06:46:23 +0000