Skip to content

Commit

Permalink
local element size for time-stepping (#279)
Browse files Browse the repository at this point in the history 
Fixes to use local element size in dt calculation

* switched to local hmin for dt calc
* removed extra 1/order in local hmin usage in dt calc
* switch to printing CFL over dt when dt is set to constant
* fixes to cfl output when constant dt is selected
* working thru pr warnings
* removed comments
* Fix maybe-uninitialized warning from gnu12

---------

Co-authored-by: Sigfried Haering <shaering@oden.utexas.edu>
Co-authored-by: Todd A. Oliver <oliver@oden.utexas.edu>
  • Loading branch information
@shaering @trevilo
3 people authored May 29, 2024
1 parent 7a68c21 commit fff2497
Show file tree
Hide file tree
Showing 3 changed files with 45 additions and 105 deletions.
2 changes: 1 addition & 1 deletion src/algebraicSubgridModels.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,7 +114,7 @@ void AlgebraicSubgridModels::initializeSelf() {
// gridScaleZ_.SetSize(sfes_truevsize);
// resolution_gf_.SetSpace(sfes_);

if (verbose) grvy_printf(ginfo, "AlgebraicSubgridModels vectors and gf initialized...\n");
if (rank0_) grvy_printf(ginfo, "AlgebraicSubgridModels vectors and gf initialized...\n");

// exports
toFlow_interface_.eddy_viscosity = &subgridVisc_gf_;
Expand Down
146 changes: 43 additions & 103 deletions src/loMach.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -202,9 +202,7 @@ void LoMachSolver::initialize() {
temporal_coeff_.dt = 1.0;
temporal_coeff_.dt3 = temporal_coeff_.dt2 = temporal_coeff_.dt1 = temporal_coeff_.dt;
SetTimeIntegrationCoefficients(0);

CFL = loMach_opts_.ts_opts_.cfl_;
if (verbose) grvy_printf(ginfo, "got CFL...\n");
CFL_ = loMach_opts_.ts_opts_.cfl_;

//-----------------------------------------------------------
// NOTE: Aspects of the ordering below are important. Beware when
Expand Down Expand Up @@ -346,12 +344,21 @@ void LoMachSolver::solveBegin() {
std::vector<double> flow_screen_values;
flow_->screenValues(flow_screen_values);

double max_cfl = -1;
if (loMach_opts_.ts_opts_.constant_dt_ > 0.0) {
max_cfl = computeCFL();
}

if (rank0_) {
// TODO(trevilo): Add state summary
std::cout << std::endl;
std::cout << std::setw(11) << "# Step ";
std::cout << std::setw(13) << "Time ";
std::cout << std::setw(13) << "dt ";
if (loMach_opts_.ts_opts_.constant_dt_ > 0.0) {
std::cout << std::setw(13) << "CFL ";
} else {
std::cout << std::setw(13) << "dt ";
}
std::cout << std::setw(13) << "Wtime/Step ";
for (size_t i = 0; i < thermo_header.size(); i++) {
std::cout << std::setw(12) << thermo_header[i] << " ";
Expand All @@ -365,7 +372,11 @@ void LoMachSolver::solveBegin() {

std::cout << std::setw(10) << iter << " ";
std::cout << std::setw(10) << std::scientific << temporal_coeff_.time << " ";
std::cout << std::setw(10) << std::scientific << temporal_coeff_.dt << " ";
if (loMach_opts_.ts_opts_.constant_dt_ > 0.0) {
std::cout << std::setw(10) << std::scientific << max_cfl << " ";
} else {
std::cout << std::setw(10) << std::scientific << temporal_coeff_.dt << " ";
}
std::cout << std::setw(10) << std::scientific << 0.0 << " ";
for (size_t i = 0; i < thermo_screen_values.size(); i++) {
std::cout << std::setw(10) << std::scientific << thermo_screen_values[i] << " ";
Expand Down Expand Up @@ -418,11 +429,20 @@ void LoMachSolver::solveStep() {
std::vector<double> flow_screen_values;
flow_->screenValues(flow_screen_values);

double max_cfl = -1;
if (loMach_opts_.ts_opts_.constant_dt_ > 0.0) {
max_cfl = computeCFL();
}

if (rank0_) {
// TODO(trevilo): Add state summary
std::cout << std::setw(10) << iter << " ";
std::cout << std::setw(10) << std::scientific << temporal_coeff_.time << " ";
std::cout << std::setw(10) << std::scientific << temporal_coeff_.dt << " ";
if (loMach_opts_.ts_opts_.constant_dt_ > 0.0) {
std::cout << std::setw(10) << std::scientific << max_cfl << " ";
} else {
std::cout << std::setw(10) << std::scientific << temporal_coeff_.dt << " ";
}
std::cout << std::setw(10) << std::scientific << max_time_per_step << " ";
for (size_t i = 0; i < thermo_screen_values.size(); i++) {
std::cout << std::setw(10) << std::scientific << thermo_screen_values[i] << " ";
Expand Down Expand Up @@ -515,11 +535,8 @@ void LoMachSolver::updateTimestep() {
// int Sdof = (turbModel_->getGridScale())->Size();
auto dataU = flow_->getCurrentVelocity()->HostRead();

// come in divided by order
// const double *dataD = bufferGridScale->HostRead();
// const double *dataD = (turbModel_->getGridScale())->HostRead();
// comes in divided by order
const double *dataD = (meshData_->getGridScale())->HostRead();
// int Sdof = dataD->Size();
int Sdof = meshData_->getDofSize();

for (int n = 0; n < Sdof; n++) {
Expand All @@ -533,41 +550,7 @@ void LoMachSolver::updateTimestep() {
}

MPI_Allreduce(&Umax_lcl, &max_speed, 1, MPI_DOUBLE, MPI_MAX, pmesh_->GetComm());
double dtInst_conv = 0.5 * CFL / max_speed;

// double *dGradU = flowClass->gradU.HostReadWrite();
// double *dGradV = flowClass->gradV.HostReadWrite();
// double *dGradW = flowClass->gradW.HostReadWrite();
// double *dVisc = tcClass->viscSml.HostReadWrite();
// double *rho = tcClass->rn.HostReadWrite();
// Umax_lcl = 1.0e-12;
// for (int n = 0; n < SdofInt; n++) {
// DenseMatrix gU;
// gU.SetSize(nvel_, dim_);
// for (int dir = 0; dir < dim_; dir++) {
// gU(0, dir) = dGradU[n + dir * SdofInt];
// }
// for (int dir = 0; dir < dim_; dir++) {
// gU(1, dir) = dGradV[n + dir * SdofInt];
// }
// for (int dir = 0; dir < dim_; dir++) {
// gU(2, dir) = dGradW[n + dir * SdofInt];
// }
// double duMag = 0.0;
// for (int dir = 0; dir < dim_; dir++) {
// for (int eq = 0; eq < nvel_; eq++) {
// duMag += gU(eq, dir) * gU(eq, dir);
// }
// }
// duMag = sqrt(duMag);
// double viscVel = sqrt(duMag * dVisc[n] / rho[n]);
// viscVel /= dataD[n];
// Umax_lcl = std::max(viscVel, Umax_lcl);
// }
// MPI_Allreduce(&Umax_lcl, &max_speed, 1, MPI_DOUBLE, MPI_MAX, pmesh_->GetComm());
// double dtInst_visc = 0.5 * CFL / max_speed;

// double dtInst = max(dtInst_conv, dtInst_visc);
double dtInst_conv = CFL_ / std::max(max_speed, 1.0e-12);

double dtInst = dtInst_conv;
double &dt = temporal_coeff_.dt;
Expand All @@ -592,21 +575,16 @@ void LoMachSolver::updateTimestep() {

double LoMachSolver::computeCFL() {
const double dt = temporal_coeff_.dt;

double Umax_lcl = 1.0e-12;
double max_speed = Umax_lcl;
double Umag;
// int Sdof = sfes_->GetNDofs();
// int Sdof = (turbModel_->getGridScale())->Size();

// come in divided by order
// double *dataD = bufferGridScale->HostReadWrite();
// double *dataD = (turbModel_->getGridScale())->HostReadWrite();
// comes in divided by order
auto dataU = flow_->getCurrentVelocity()->HostRead();
const double *dataD = (meshData_->getGridScale())->HostRead();
// int Sdof = dataD->Size();
int Sdof = meshData_->getDofSize();

MPI_Barrier(groupsMPI->getTPSCommWorld());
for (int n = 0; n < Sdof; n++) {
Umag = 0.0;
// Vector delta({dataX[n], dataY[n], dataZ[n]});
Expand All @@ -617,60 +595,20 @@ double LoMachSolver::computeCFL() {
Umag = std::sqrt(Umag);
Umax_lcl = std::max(Umag, Umax_lcl);
}
MPI_Allreduce(&Umax_lcl, &max_speed, 1, MPI_DOUBLE, MPI_MAX, pmesh_->GetComm());
double CFL_conv = 2.0 * dt * max_speed;

// double *dGradU = flowClass->gradU.HostReadWrite();
// double *dGradV = flowClass->gradV.HostReadWrite();
// double *dGradW = flowClass->gradW.HostReadWrite();
// double *dVisc = tcClass->viscSml.HostReadWrite();
// double *rho = tcClass->rn.HostReadWrite();
// Umax_lcl = 1.0e-12;
// for (int n = 0; n < SdofInt; n++) {
// DenseMatrix gU;
// gU.SetSize(nvel_, dim_);
// for (int dir = 0; dir < dim_; dir++) {
// gU(0, dir) = dGradU[n + dir * SdofInt];
// }
// for (int dir = 0; dir < dim_; dir++) {
// gU(1, dir) = dGradV[n + dir * SdofInt];
// }
// for (int dir = 0; dir < dim_; dir++) {
// gU(2, dir) = dGradW[n + dir * SdofInt];
// }
// double duMag = 0.0;
// for (int dir = 0; dir < dim_; dir++) {
// for (int eq = 0; eq < nvel_; eq++) {
// duMag += gU(eq, dir) * gU(eq, dir);
// }
// }
// duMag = sqrt(duMag);
// double viscVel = sqrt(duMag * dVisc[n] / rho[n]);
// viscVel /= dataD[n];
// Umax_lcl = std::max(viscVel, Umax_lcl);
// }
// MPI_Allreduce(&Umax_lcl, &max_speed, 1, MPI_DOUBLE, MPI_MAX, pmesh_->GetComm());
// double CFL_visc = 2.0 * dt * max_speed;

// double CFL_here = max(CFL_conv, CFL_visc);
double CFL_here = CFL_conv;

return CFL_here;
MPI_Allreduce(&Umax_lcl, &max_speed, 1, MPI_DOUBLE, MPI_MAX, groupsMPI->getTPSCommWorld());
double CFL_conv = dt * max_speed;

return CFL_conv;
}

void LoMachSolver::setTimestep() {
double Umax_lcl = 1.0e-12;
double max_speed = Umax_lcl;
double convT_lcl = 1.0e-12;
double min_convT = 1.0;
double Umag;
// int Sdof = sfes_->GetNDofs();
// int Sdof = (turbModel_->getGridScale())->Size();
// const double *dataD = (meshData_->getGridScale())->HostRead();
// int Sdof = dataD->Size();
hmin = meshData_->getMinGridScale();
const double *dataD = (meshData_->getGridScale())->HostRead();
int Sdof = meshData_->getDofSize();

CFL = loMach_opts_.ts_opts_.cfl_;

// dt_fixed is initialized to -1, so if it is positive,
// then the user requested a fixed dt run
const double dt_fixed = loMach_opts_.ts_opts_.constant_dt_;
Expand All @@ -686,17 +624,19 @@ void LoMachSolver::setTimestep() {
}
Umag = std::sqrt(Umag);
Umax_lcl = std::max(Umag, Umax_lcl);

// local hmin comes in devided by order
double hmin = dataD[n];
convT_lcl = hmin / Umax_lcl;
}
MPI_Allreduce(&Umax_lcl, &max_speed, 1, MPI_DOUBLE, MPI_MAX, pmesh_->GetComm());
double dtInst = CFL * hmin / (max_speed * (double)order);
MPI_Allreduce(&convT_lcl, &min_convT, 1, MPI_DOUBLE, MPI_MIN, pmesh_->GetComm());
double dtInst = CFL_ * min_convT;
temporal_coeff_.dt = dtInst;

const double dt_initial = loMach_opts_.ts_opts_.initial_dt_;
if ((dt_initial < dtInst) && !(loMach_opts_.io_opts_.enable_restart_)) {
temporal_coeff_.dt = dt_initial;
}

std::cout << "dt from setTimestep: " << temporal_coeff_.dt << " max_speed: " << max_speed << endl;
}
}

Expand Down
2 changes: 1 addition & 1 deletion src/loMach.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -180,7 +180,7 @@ class LoMachSolver : public TPS::PlasmaSolver {
// Time marching related parameters
int iter;
int iter_start_;
double CFL;
double CFL_;

int max_bdf_order; // input option now
temporalSchemeCoefficients temporal_coeff_;
Expand Down

0 comments on commit fff2497

Please sign in to comment.