diff options
Diffstat (limited to 'src/nullsurf.c')
| -rw-r--r-- | src/nullsurf.c | 38 |
1 files changed, 26 insertions, 12 deletions
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; } |