From 1253302aaf002925c76e5ab0f58c4c5c961a2e54 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Sat, 30 May 2020 12:05:51 +0200 Subject: Implement multiple surfaces. --- interface.ccl | 7 ++++--- param.ccl | 5 +++++ schedule.ccl | 3 ++- src/nullsurf.c | 38 ++++++++++++++++++++++++++------------ 4 files changed, 37 insertions(+), 16 deletions(-) diff --git a/interface.ccl b/interface.ccl index b0f1325..18580c9 100644 --- a/interface.ccl +++ b/interface.ccl @@ -15,7 +15,7 @@ USES FUNCTION VarDataPtrI public: -REAL photon_coord TYPE=array DIM=2 SIZE=nb_surfaces,photons_per_surface DISTRIB=constant +REAL photon_coord TYPE=array DIM=2 SIZE=photons_per_surface,nb_surfaces DISTRIB=constant { pos_x, # stored only to send zeroes to interpolator @@ -25,7 +25,7 @@ REAL photon_coord TYPE=array DIM=2 SIZE=nb_surfaces,photons_per_surface DISTRIB= mom_z } -REAL photon_rhs TYPE=array DIM=3 SIZE=4,nb_surfaces,photons_per_surface DISTRIB=constant +REAL photon_rhs TYPE=array DIM=3 SIZE=photons_per_surface,nb_surfaces,4 DISTRIB=constant { pos_x_rhs, pos_z_rhs, @@ -33,6 +33,7 @@ REAL photon_rhs TYPE=array DIM=3 SIZE=4,nb_surfaces,photons_per_surface DISTRIB= mom_z_rhs } -REAL photon_times TYPE=array DIM=1 SIZE=nb_surfaces DISTRIB=constant INT evol_next_step TYPE=scalar INT surfaces_active TYPE=scalar +REAL origin_proper_time TYPE=scalar +REAL origin_proper_time_rhs TYPE=array DIM=1 SIZE=4 diff --git a/param.ccl b/param.ccl index af9189c..25286ad 100644 --- a/param.ccl +++ b/param.ccl @@ -10,6 +10,11 @@ CCTK_INT nb_surfaces "Number of null surfaces to trace" 1: :: "" } 1 +CCTK_REAL launch_delta "Proper time between subsequent surface launches" +{ + 0: :: "" +} 0.25 + CCTK_INT solve_level "Refinement level to integrate photons on" { 0: :: "" diff --git a/schedule.ccl b/schedule.ccl index e743105..9b32340 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -16,7 +16,8 @@ if (nb_surfaces > 0) { STORAGE: photon_coord STORAGE: photon_rhs - STORAGE: photon_times STORAGE: evol_next_step STORAGE: surfaces_active + STORAGE: origin_proper_time + STORAGE: origin_proper_time_rhs } diff --git a/src/nullsurf.c b/src/nullsurf.c index 11563ec..b97ee44 100644 --- a/src/nullsurf.c +++ b/src/nullsurf.c @@ -473,8 +473,16 @@ ns_evol_rk_substep(NSContext *ns, int substep, int surfaces_active) double *mom_x = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::mom_x"); double *mom_z = CCTK_VarDataPtr(ns->gh, 0, "NullSurf::mom_z"); + const int origin_idx = CCTK_GFINDEX3D(ns->gh, ns->gh->cctk_nghostzones[0], + 0, ns->gh->cctk_nghostzones[2]); + double *tau_rhs = (double*)CCTK_VarDataPtr(ns->gh, 0, "NullSurf::origin_proper_time_rhs") + substep; + double alpha = ((double*)CCTK_VarDataPtr(ns->gh, 0, metric_vars[ALPHA]))[origin_idx]; + double *step_pos_x, *step_pos_z, *step_mom_x, *step_mom_z; + /* proper time at origin */ + *tau_rhs = alpha; + /* compute the coordinate/momenta to evaluate the RHS for */ switch (substep) { case 0: @@ -548,6 +556,11 @@ static void ns_evol_rk_finish(NSContext *ns, int surfaces_active) const double *mom_z_rhs2 = mom_z_rhs + 2 * rk_rhs_stride; const double *mom_z_rhs3 = mom_z_rhs + 3 * rk_rhs_stride; + double *tau_rhs = (double*)CCTK_VarDataPtr(ns->gh, 0, "NullSurf::origin_proper_time_rhs"); + double *tau = (double*)CCTK_VarDataPtr(ns->gh, 0, "NullSurf::origin_proper_time"); + + *tau += ns->dt * (tau_rhs[0] / 6. + tau_rhs[1] / 3. + tau_rhs[2] / 3. + tau_rhs[3] / 6.); + #pragma omp parallel for for (int i = 0; i < surfaces_active * ns->photons_per_surface; i++) { pos_x[i] += ns->dt * (pos_x_rhs0[i] / 6. + pos_x_rhs1[i] / 3. + pos_x_rhs2[i] / 3. + pos_x_rhs3[i] / 6.); @@ -573,24 +586,28 @@ void ns_evol(CCTK_ARGUMENTS) fprintf(stderr, "ns evol step %d at time %g; pos[zaxis] = %g; pos[xaxis] = %g\n", *evol_next_step, cctk_time, pos_z[0], pos_x[photons_per_surface - 1]); + /* do the last Runge-Kutta sub-step */ + if (*evol_next_step == 3) { + ns_evol_rk_substep(ns, *evol_next_step, *surfaces_active); + ns_evol_rk_finish(ns, *surfaces_active); + *evol_next_step = 0; + } + /* launch a new surface if its time has come */ - // TODO: actually implement multiple surfaces - if (cctk_time == 0.0 && *surfaces_active == 0) { + if (*surfaces_active < nb_surfaces && + *origin_proper_time >= *surfaces_active * launch_delta) { const int offset = *surfaces_active * photons_per_surface; + fprintf(stderr, "launching surface %d at t=%g/tau=%g\n", + *surfaces_active, cctk_time, *origin_proper_time); surface_launch(ns, pos_x + offset, pos_z + offset, mom_x + offset, mom_z + offset); - *surfaces_active += 1; } - /* do the Runge-Kutta evolution step(s) */ + /* do the Runge-Kutta evolution step(s): either 0->1 or 1->2 */ ns_evol_rk_substep(ns, *evol_next_step, *surfaces_active); *evol_next_step += 1; - if (*evol_next_step == 4) { - ns_evol_rk_finish(ns, *surfaces_active); - *evol_next_step = 0; - } - if (*evol_next_step == 0 || *evol_next_step == 2) { + if (*evol_next_step == 2) { ns_evol_rk_substep(ns, *evol_next_step, *surfaces_active); *evol_next_step += 1; } @@ -686,9 +703,6 @@ void ns_init(CCTK_ARGUMENTS) if (reflevel != solve_level) return; - for (int i = 0; i < nb_surfaces; i++) - photon_times[i] = -1; - *surfaces_active = 0; *evol_next_step = 0; } -- cgit v1.2.3