Skip to content

Commit

Permalink
comments and other minor.
Browse files Browse the repository at this point in the history
  • Loading branch information
@vladotomov
vladotomov committed Dec 18, 2024
1 parent e7a7310 commit 1b148ee
Show file tree
Hide file tree
Showing 3 changed files with 49 additions and 69 deletions.
22 changes: 11 additions & 11 deletions remhos.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -333,7 +333,7 @@ int main(int argc, char *argv[])
if (ode_solver_type > 10)
{
auto rk_idp = dynamic_cast<RKIDPSolver *>(idp_ode_solver);
if (rk_idp) { rk_idp->UseMask(true); }
if (rk_idp) { rk_idp->UseMask(false); }
ode_solver = idp_ode_solver;
}

Expand Down Expand Up @@ -934,8 +934,7 @@ int main(int argc, char *argv[])

double u_min, u_max;
GetMinMax(u, u_min, u_max);
double s_min = numeric_limits<double>::infinity(),
s_max = -numeric_limits<double>::infinity();
double s_min = -1.0, s_max = -1.0;
if (product_sync) { ComputeMinMaxS(NE, us, u, s_min, s_max); }

if (exec_mode == 1)
Expand Down Expand Up @@ -1025,6 +1024,7 @@ int main(int argc, char *argv[])
// Monotonicity check for debug purposes mainly.
if (verify_bounds && forced_bounds && smth_indicator == NULL)
{
const double eps = 1e-10;
double u_min_new, u_max_new,
s_min_new = s_min, s_max_new = s_max;
GetMinMax(u, u_min_new, u_max_new);
Expand All @@ -1037,13 +1037,13 @@ int main(int argc, char *argv[])
{
if (myid == 0)
{
MFEM_VERIFY(u_min_new > u_min - 1e-12,
MFEM_VERIFY(u_min_new > u_min - eps,
"Undershoot of " << u_min - u_min_new);
MFEM_VERIFY(u_max_new < u_max + 1e-12,
MFEM_VERIFY(u_max_new < u_max + eps,
"Overshoot of " << u_max_new - u_max);
MFEM_VERIFY(s_min_new > s_min - 1e-12,
MFEM_VERIFY(s_min_new > s_min - eps,
"Undershoot in s of " << s_min - s_min_new);
MFEM_VERIFY(s_max_new < s_max + 1e-12,
MFEM_VERIFY(s_max_new < s_max + eps,
"Overshoot in s of " << s_max_new - s_max);
}
u_min = u_min_new;
Expand All @@ -1053,13 +1053,13 @@ int main(int argc, char *argv[])
{
if (myid == 0)
{
MFEM_VERIFY(u_min_new > 0.0 - 1e-12,
MFEM_VERIFY(u_min_new > 0.0 - eps,
"Undershoot of " << 0.0 - u_min_new);
MFEM_VERIFY(u_max_new < 1.0 + 1e-12,
MFEM_VERIFY(u_max_new < 1.0 + eps,
"Overshoot of " << u_max_new - 1.0);
MFEM_VERIFY(s_min_new > 0.0 - 1e-12,
MFEM_VERIFY(s_min_new > 0.0 - eps,
"Undershoot in s of " << 0.0 - s_min_new);
MFEM_VERIFY(s_max_new < 1.0 + 1e-12,
MFEM_VERIFY(s_max_new < 1.0 + eps,
"Overshoot in s of " << s_max_new - 1.0);
}
}
Expand Down
89 changes: 34 additions & 55 deletions remhos_fct.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,32 +115,6 @@ void FCTSolver::CalcCompatibleLOProduct(const ParGridFunction &us,
dof_id = k*ndofs + j;
d_us_LO_new(dof_id) = (u_new(dof_id) * s_avg - us(dof_id)) / dt;
}

#ifdef REMHOS_FCT_PRODUCT_DEBUG
// Check the LO product solution.
double us_min, us_max;
for (int j = 0; j < ndofs; j++)
{
dof_id = k*ndofs + j;
if (active_dofs[dof_id] == false) { continue; }

us_min = s_min_loc(j) * u_new(dof_id);
us_max = s_max_loc(j) * u_new(dof_id);

if (s_avg * u_new(dof_id) + eps < us_min ||
s_avg * u_new(dof_id) - eps > us_max)
{
std::cout << "---\ns_avg * u: " << k << " "
<< us_min << " "
<< s_avg * u_new(dof_id) << " "
<< us_max << endl
<< u_new(dof_id) << " " << s_avg << endl
<< s_min_loc(j) << " " << s_max_loc(j) << "\n---\n";

MFEM_ABORT("s_avg * u not in bounds");
}
}
#endif
}
}

Expand Down Expand Up @@ -287,43 +261,48 @@ void FluxBasedFCT::CalcFCTProduct(const ParGridFunction &us, const Vector &m,
dus_lo_fct = d_us;
}

#ifdef REMHOS_FCT_PRODUCT_DEBUG
// Check the bounds of the final solution.
const double eps = 1e-12;
Vector us_new(d_us.Size());
add(1.0, us, dt, d_us, us_new);
for (int k = 0; k < NE; k++)
if (verify_bounds)
{
if (active_el[k] == false) { continue; }

for (int j = 0; j < ndofs; j++)
// Check the bounds of the final solution.
const double eps = 1e-12;
Vector us_new(d_us.Size());
add(1.0, us, dt, d_us, us_new);
for (int k = 0; k < NE; k++)
{
dof_id = k*ndofs + j;
if (active_dofs[dof_id] == false) { continue; }
if (active_el[k] == false) { continue; }

double s = us_new(dof_id) / u_new(dof_id);
if (s + eps < s_min(dof_id) ||
s - eps > s_max(dof_id))
for (int j = 0; j < ndofs; j++)
{
std::cout << "Final s " << j << " " << k << " "
<< s_min(dof_id) << " "
<< s << " "
<< s_max(dof_id) << std::endl;
std::cout << "---\n";
}
dof_id = k*ndofs + j;
if (active_dofs[dof_id] == false) { continue; }

if (us_new(dof_id) + eps < us_min(dof_id) ||
us_new(dof_id) - eps > us_max(dof_id))
{
std::cout << "Final us " << j << " " << k << " "
<< us_min(dof_id) << " "
<< us_new(dof_id) << " "
<< us_max(dof_id) << std::endl;
std::cout << "---\n";
double s = us_new(dof_id) / u_new(dof_id);
if (s + eps < s_min(dof_id) ||
s - eps > s_max(dof_id))
{
std::cout << "Final s " << j << " " << k << " "
<< s_min(dof_id) << " "
<< s << " "
<< s_max(dof_id) << std::endl;
std::cout << "---\n";

MFEM_ABORT("s not in bounds after FCT.");
}

if (us_new(dof_id) + eps < us_min(dof_id) ||
us_new(dof_id) - eps > us_max(dof_id))
{
std::cout << "Final us " << j << " " << k << " "
<< us_min(dof_id) << " "
<< us_new(dof_id) << " "
<< us_max(dof_id) << std::endl;
std::cout << "---\n";

MFEM_ABORT("us not in bounds after FCT.");
}
}
}
}
#endif
}

void FluxBasedFCT::ComputeFluxMatrix(const ParGridFunction &u,
Expand Down
7 changes: 4 additions & 3 deletions remhos_solvers.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -212,25 +212,26 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt)
// Update mask with the HO update
if (use_masks) { UpdateMask(x, dxs[i], dct, mask); }

//
// Form the unlimited update for the stage.
// Note that it converts eq. (2.16) in JLG's paper into an update using
// the previous limited updates.
if (d_i[i] != 1.)
//
{
// for mask = 0, we get dxs (nothing happens).
// the loop below won't change it -> Forward Euler.
// for mask = 1, we scale dxs by d_i[i].
if (use_masks) { AddMasked(mask, d_i[i]-1., dxs[i], dxs[i]); }
else { dxs[i] *= d_i[i]; }
}
// Use all previous limited updates.
for (int j = 0; j < i; j++)
{
// Use all previous limited updates.
if (use_masks) { AddMasked(mask, d_i[j], dxs[j], dxs[i]); }
else { dxs[i].Add(d_i[j], dxs[j]); }
}

// Limit the step
// Limit the step (always a Forward Euler step).
f->LimitMult(x, dxs[i]);

// Update the state
Expand Down

0 comments on commit 1b148ee

Please sign in to comment.