Skip to content

Commit

Permalink
timestepping problem is somewhat solved. In the future, there should …
Browse files Browse the repository at this point in the history 
…be per-block dt_last instead.
  • Loading branch information
Hyerin Cho committed Jun 24, 2024
1 parent 1bbb2eb commit a5fad8b
Show file tree
Hide file tree
Showing 4 changed files with 221 additions and 42 deletions.
6 changes: 3 additions & 3 deletions kharma/driver/kharma_driver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -469,9 +469,9 @@ TaskID KHARMADriver::AddStateUpdateIdealGuess(TaskID& t_start, TaskList& tl, Mes
void KHARMADriver::SetGlobalTimeStep()
{
// TODO(BSP) apply the limits from GRMHD package here
if (tm.dt < 0.1 * std::numeric_limits<Real>::max()) {
tm.dt *= 2.0;
}
//if (tm.dt < 0.1 * std::numeric_limits<Real>::max()) {
// tm.dt *= 2.0;
//}
Real big = std::numeric_limits<Real>::max();
for (auto const &pmb : pmesh->block_list) {
tm.dt = std::min(tm.dt, pmb->NewDt());
Expand Down
49 changes: 42 additions & 7 deletions kharma/driver/multizone_step.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,30 +68,33 @@ TaskCollection KHARMADriver::MakeMultizoneTaskCollection(BlockList_t &blocks, in
bool is_active[num_blocks]; // = {false, false, true};
bool apply_boundary_condition[num_blocks][BOUNDARY_NFACES];

Multizone::DecideActiveBlocks(pmesh, tm, is_active, apply_boundary_condition);
Multizone::DecideActiveBlocksAndBoundaryConditions(pmesh, tm, is_active, apply_boundary_condition);
for (int i = 0; i < num_partitions; i++)
std::cout << "iblock " << i << ": is active? " << is_active[i] << ", boundary applied? " << apply_boundary_condition[i][BoundaryFace::inner_x1] << apply_boundary_condition[i][BoundaryFace::outer_x1] << std::endl;

// TaskCollections are a collection of TaskRegions.
// Each TaskRegion can operate on eash meshblock separately, i.e. one MeshBlockData object (slower),
// or on a collection of MeshBlock objects called the MeshData
TaskCollection tc;
const TaskID t_none(0);

Flag("MakeTaskCollection::timestep");
//Flag("MakeTaskCollection::timestep");

// Timestep region: calculate timestep based on the newly updated active zones
TaskRegion &timestep_region = tc.AddRegion(num_partitions);
//TaskRegion &timestep_region = tc.AddRegion(num_partitions);
// Estimate next time step based on ctop
for (int i = 0; i < num_partitions; i++) {
auto &tl = timestep_region[i];
//auto &tl = timestep_region[i];
if (is_active[i]) {
auto &base = pmesh->mesh_data.GetOrAdd("base", i);
auto t_new_dt =
tl.AddTask(t_none, Update::EstimateTimestep<MeshData<Real>>, base.get());
Update::EstimateTimestep<MeshData<Real>>(base.get());
//auto t_new_dt =
// tl.AddTask(t_none, Update::EstimateTimestep<MeshData<Real>>, base.get());
}
}
SetGlobalTimeStep();

EndFlag();
//EndFlag();

// Which packages we load affects which tasks we'll add to the list
auto& pkgs = pmesh->packages.AllPackages();
Expand Down Expand Up @@ -302,6 +305,38 @@ TaskCollection KHARMADriver::MakeMultizoneTaskCollection(BlockList_t &blocks, in
KHARMADriver::AddFullSyncRegion(tc, md_sync);
}

EndFlag();

Flag("MakeTaskCollection::timestep");
// HYERIN (06/20/24) splitting this part to the end for now
// Switch region: decide if we want to switch zones
TaskRegion &switch_region = tc.AddRegion(1);
auto &tl = switch_region[0];
auto &md_temp = pmesh->mesh_data.Add(integrator->stage_name[stage]);
bool switch_zone = false;
auto t_switch = t_none;
if (integrator->nstages == stage) {
t_switch = tl.AddTask(t_none, Multizone::DecideToSwitch, md_temp.get(), tm, switch_zone);
}

// Timestep region: calculate timestep based on the newly updated active zones
//TaskRegion &timestep_region = tc.AddRegion(num_partitions);
//auto t_new_active = t_none;
//// Estimate next time step based on ctop
//for (int i = 0; i < num_partitions; i++) {
// auto &tl = timestep_region[i];
// //if (switch_zone) { // take a next step and re-evaluate the active block
// t_new_active = tl.AddTask(t_none, Multizone::DecideNextActiveBlocks, md_temp.get(), tm, i, is_active[i], switch_zone);
// //}
// if (is_active[i]) {
// auto &base = pmesh->mesh_data.GetOrAdd("base", i);
// // Update::EstimateTimestep<MeshData<Real>>(base.get());
// auto t_new_dt =
// tl.AddTask(t_new_active, Update::EstimateTimestep<MeshData<Real>>, base.get());
// }
//}


EndFlag();

return tc;
Expand Down
191 changes: 162 additions & 29 deletions kharma/multizone/multizone.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,11 +80,11 @@ std::shared_ptr<KHARMAPackage> Multizone::Initialize(ParameterInput *pin, std::s
return pkg;
}

void Multizone::DecideActiveBlocks(Mesh *pmesh, const SimTime &tm, bool *is_active, bool apply_boundary_condition[][BOUNDARY_NFACES])
void Multizone::DecideActiveBlocksAndBoundaryConditions(Mesh *pmesh, const SimTime &tm, bool *is_active, bool apply_boundary_condition[][BOUNDARY_NFACES])
{
Flag("DecideActiveBlocks");
Flag("DecideActiveBlocksAndBoundaryConditions");
auto &params = pmesh->packages.Get("Multizone")->AllParams();

// Input parameters
const int nzones = params.Get<int>("nzones");
const Real base = params.Get<Real>("base");
Expand All @@ -99,33 +99,33 @@ void Multizone::DecideActiveBlocks(Mesh *pmesh, const SimTime &tm, bool *is_acti
// Current location in V-cycles
int i_vcycle = params.Get<int>("i_vcycle"); // current i_vcycle
int i_within_vcycle = params.Get<int>("i_within_vcycle"); // current i_within_vcycle
Real t0_zone = params.Get<Real>("t0_zone"); // time at zone-switching
int n0_zone = params.Get<int>("n0_zone"); // cycle at zone-switching
//Real t0_zone = params.Get<Real>("t0_zone"); // time at zone-switching
//int n0_zone = params.Get<int>("n0_zone"); // cycle at zone-switching
int i_zone = m::abs(i_within_vcycle - (nzones - 1)); // zone number. zone-0 is the smallest annulus.

// Derived parameters
bool switch_zone = false; // determines if we are switching zones here
Real temp_rin, runtime_per_zone;
int nzones_per_vcycle = 2 * (nzones - 1);
if (ncycle_per_zone > 0) switch_zone = (tm.ncycle - n0_zone) >= ncycle_per_zone;
else {
temp_rin = m::pow(base, i_zone);
runtime_per_zone = f_tchar * CalcRuntime(temp_rin, base, gam, bondi_rs, loc_tchar);
switch_zone = (tm.time - t0_zone) > runtime_per_zone;
//std::cout << "time now " << tm.time << ", t0_zone " << t0_zone << ", runtime_per_zone " << runtime_per_zone << std::endl;
}
if (switch_zone) {
i_within_vcycle += 1;
if (i_within_vcycle >= nzones_per_vcycle) { // if completed a V-cycle
i_within_vcycle -= nzones_per_vcycle;
i_vcycle += 1;
}
i_zone = m::abs(i_within_vcycle - (nzones - 1)); // update zone number
params.Update<int>("i_within_vcycle", i_within_vcycle);
params.Update<int>("i_vcycle", i_vcycle);
params.Update<Real>("t0_zone", tm.time);
params.Update<int>("n0_zone", tm.ncycle);
}
// Derived parameters (split into other function)
//bool switch_zone = false; // determines if we are switching zones here
//Real temp_rin, runtime_per_zone;
//int nzones_per_vcycle = 2 * (nzones - 1);
//if (ncycle_per_zone > 0) switch_zone = (tm.ncycle - n0_zone) >= ncycle_per_zone;
//else {
// temp_rin = m::pow(base, i_zone);
// runtime_per_zone = f_tchar * CalcRuntime(temp_rin, base, gam, bondi_rs, loc_tchar);
// switch_zone = (tm.time - t0_zone) > runtime_per_zone;
// //std::cout << "time now " << tm.time << ", t0_zone " << t0_zone << ", runtime_per_zone " << runtime_per_zone << std::endl;
//}
//if (switch_zone) {
// i_within_vcycle += 1;
// if (i_within_vcycle >= nzones_per_vcycle) { // if completed a V-cycle
// i_within_vcycle -= nzones_per_vcycle;
// i_vcycle += 1;
// }
// i_zone = m::abs(i_within_vcycle - (nzones - 1)); // update zone number
// params.Update<int>("i_within_vcycle", i_within_vcycle);
// params.Update<int>("i_vcycle", i_vcycle);
// params.Update<Real>("t0_zone", tm.time);
// params.Update<int>("n0_zone", tm.ncycle);
//}

// Range of radii that is active
int active_rout;
Expand Down Expand Up @@ -157,10 +157,143 @@ void Multizone::DecideActiveBlocks(Mesh *pmesh, const SimTime &tm, bool *is_acti
if (m::abs(x1min - active_x1min) / m::max(active_x1min, SMALL) < 1.e-10) apply_boundary_condition[iblock][BoundaryFace::inner_x1] = true;
if (m::abs(x1max - active_x1max) / active_x1max < 1.e-10) apply_boundary_condition[iblock][BoundaryFace::outer_x1] = true;
}
std::cout << "iblock " << iblock << " x1min " << x1min << " x1max " << x1max << ": is active? " << is_active[iblock] << ", boundary applied? " << apply_boundary_condition[iblock][BoundaryFace::inner_x1] << apply_boundary_condition[iblock][BoundaryFace::outer_x1] << std::endl;
//std::cout << "iblock " << iblock << " x1min " << x1min << " x1max " << x1max << ": is active? " << is_active[iblock] << ", boundary applied? " << apply_boundary_condition[iblock][BoundaryFace::inner_x1] << apply_boundary_condition[iblock][BoundaryFace::outer_x1] << std::endl;
}
EndFlag();

}

TaskStatus Multizone::DecideToSwitch(MeshData<Real> *md, const SimTime &tm, bool &switch_zone)
{
Flag("DecideToSwitch");
auto pmesh = md->GetMeshPointer();
auto &params = pmesh->packages.Get("Multizone")->AllParams();

const int next_ncycle = tm.ncycle + 1;
const Real next_time = tm.time + tm.dt;

// Input parameters
const int nzones = params.Get<int>("nzones");
const Real base = params.Get<Real>("base");
const int num_Vcycles = params.Get<int>("num_Vcycles");
const int ncycle_per_zone = params.Get<int>("ncycle_per_zone");
const Real bondi_rs = params.Get<Real>("bondi_rs");
const Real gam = pmesh->packages.Get("GRMHD")->Param<Real>("gamma");
const Real f_tchar = params.Get<Real>("f_tchar");
const bool loc_tchar = params.Get<bool>("loc_tchar");

// Current location in V-cycles
int i_vcycle = params.Get<int>("i_vcycle"); // current i_vcycle
int i_within_vcycle = params.Get<int>("i_within_vcycle"); // current i_within_vcycle
Real t0_zone = params.Get<Real>("t0_zone"); // time at zone-switching
int n0_zone = params.Get<int>("n0_zone"); // cycle at zone-switching
int i_zone = m::abs(i_within_vcycle - (nzones - 1)); // zone number. zone-0 is the smallest annulus.

// Derived parameters
switch_zone = false; // default
Real temp_rin, runtime_per_zone;
int nzones_per_vcycle = 2 * (nzones - 1);
if (ncycle_per_zone > 0) switch_zone = (next_ncycle - n0_zone) >= ncycle_per_zone;
else {
temp_rin = m::pow(base, i_zone);
runtime_per_zone = f_tchar * CalcRuntime(temp_rin, base, gam, bondi_rs, loc_tchar);
switch_zone = (next_time - t0_zone) > runtime_per_zone;
//std::cout << "time now " << tm.time << ", t0_zone " << t0_zone << ", runtime_per_zone " << runtime_per_zone << std::endl;
}
if (switch_zone) {
i_within_vcycle += 1;
if (i_within_vcycle >= nzones_per_vcycle) { // if completed a V-cycle
i_within_vcycle -= nzones_per_vcycle;
i_vcycle += 1;
}
i_zone = m::abs(i_within_vcycle - (nzones - 1)); // update zone number
params.Update<int>("i_within_vcycle", i_within_vcycle);
params.Update<int>("i_vcycle", i_vcycle);
params.Update<Real>("t0_zone", next_time);
params.Update<int>("n0_zone", next_ncycle);
printf("HYERIN2: switching zone next!\n");
}

EndFlag();
return TaskStatus::complete;
}

TaskStatus Multizone::DecideNextActiveBlocks(MeshData<Real> *md, const SimTime &tm, const int iblock, bool &is_active, const bool switch_zone)
{
Flag("DecideNextActiveBlocks");
if (switch_zone) {
printf("HYERIN: switching zone here! the block is for %d \n", iblock);
auto pmesh = md->GetMeshPointer();
auto &params = pmesh->packages.Get("Multizone")->AllParams();

const int next_ncycle = tm.ncycle + 1;
const Real next_time = tm.time + tm.dt;

// Input parameters
const int nzones = params.Get<int>("nzones");
const Real base = params.Get<Real>("base");
const int num_Vcycles = params.Get<int>("num_Vcycles");
//const int ncycle_per_zone = params.Get<int>("ncycle_per_zone");
const bool move_rin = params.Get<bool>("move_rin");
//const Real bondi_rs = params.Get<Real>("bondi_rs");
//const Real gam = pmesh->packages.Get("GRMHD")->Param<Real>("gamma");
//const Real f_tchar = params.Get<Real>("f_tchar");
//const bool loc_tchar = params.Get<bool>("loc_tchar");

// Current location in V-cycles
//int i_vcycle = params.Get<int>("i_vcycle"); // current i_vcycle
int i_within_vcycle = params.Get<int>("i_within_vcycle"); // current i_within_vcycle
//Real t0_zone = params.Get<Real>("t0_zone"); // time at zone-switching
//int n0_zone = params.Get<int>("n0_zone"); // cycle at zone-switching
int i_zone = m::abs(i_within_vcycle - (nzones - 1)); // zone number. zone-0 is the smallest annulus.

// Derived parameters
//bool switch_zone = false; // determines if we are switching zones here
//Real temp_rin, runtime_per_zone;
//int nzones_per_vcycle = 2 * (nzones - 1);
//if (ncycle_per_zone > 0) switch_zone = (next_ncycle - n0_zone) >= ncycle_per_zone;
//else {
// temp_rin = m::pow(base, i_zone);
// runtime_per_zone = f_tchar * CalcRuntime(temp_rin, base, gam, bondi_rs, loc_tchar);
// switch_zone = (next_time - t0_zone) > runtime_per_zone;
// //std::cout << "time now " << tm.time << ", t0_zone " << t0_zone << ", runtime_per_zone " << runtime_per_zone << std::endl;
//}
//if (switch_zone) {
// i_within_vcycle += 1;
// if (i_within_vcycle >= nzones_per_vcycle) { // if completed a V-cycle
// i_within_vcycle -= nzones_per_vcycle;
// i_vcycle += 1;
// }
// i_zone = m::abs(i_within_vcycle - (nzones - 1)); // update zone number
// params.Update<int>("i_within_vcycle", i_within_vcycle);
// params.Update<int>("i_vcycle", i_vcycle);
// params.Update<Real>("t0_zone", next_time);
// params.Update<int>("n0_zone", next_ncycle);
//}

// Range of radii that is active
int active_rout;
int active_rin = m::pow(base, i_zone);
if (move_rin) active_rout = m::pow(base, nzones + 1);
else active_rout = m::pow(base, i_zone + 2);
Real active_x1min = m::log(active_rin);
Real active_x1max = m::log(active_rout);

const int num_blocks = pmesh->block_list.size();

// Determine Active Blocks
GReal x1min, x1max;
//for (int iblock=0; iblock < num_blocks; iblock++) {
auto &pmb = pmesh->block_list[iblock];
x1min = pmb->block_size.xmin(X1DIR);
x1max = pmb->block_size.xmax(X1DIR);
if ((x1min + 1.0e-8 < active_x1min) || (x1max - 1.0e-8 > active_x1max)) is_active = false; //[iblock]
else is_active = true;
//}
}
EndFlag();

return TaskStatus::complete;
}

TaskStatus Multizone::AverageEMFSeams(MeshData<Real> *md_emf_only, bool *apply_boundary_condition)
Expand Down
17 changes: 14 additions & 3 deletions kharma/multizone/multizone.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,11 +78,22 @@ KOKKOS_INLINE_FUNCTION Real CalcRuntime(const Real &r_in, const int &base, const
}

/**
* Decide which blocks are active
* Depending on the location in the V cycles, decide which blocks are active and where to apply boundary conditions for the current step
*
* Function in this package: Currently nothing
*/
void DecideActiveBlocks(Mesh *pmesh, const SimTime &tm, bool *is_active, bool apply_boundary_condition[][BOUNDARY_NFACES]);
void DecideActiveBlocksAndBoundaryConditions(Mesh *pmesh, const SimTime &tm, bool *is_active, bool apply_boundary_condition[][BOUNDARY_NFACES]);

/**
* Decide whether or not to progress in the V cycle
*
*/
TaskStatus DecideToSwitch(MeshData<Real> *md, const SimTime &tm, bool &switch_zone);

/**
* Decide which blocks are active for the next step and progress in the V cycle, in order to determine dt
*
*/
TaskStatus DecideNextActiveBlocks(MeshData<Real> *md, const SimTime &tm, const int iblock, bool &is_active, const bool switch_zone);

/**
* Average EMFs on seams (internal boundaries) between this block and specified other blocks.
Expand Down

0 comments on commit a5fad8b

Please sign in to comment.