Skip to content

Commit

Permalink
Semi-correct implicit solver (#21)
Browse files Browse the repository at this point in the history
This version can integrate correctly for ~200 days. But there is a memory leak now. Also, the implicit method is not exactly correct, a -> PaP^{-1}. However, it is stable without this fix.
  • Loading branch information
chengcli committed Mar 9, 2024
1 parent 78d31ee commit 1d2a847
Show file tree
Hide file tree
Showing 6 changed files with 160 additions and 31 deletions.
20 changes: 10 additions & 10 deletions examples/2023-Chen-exo3/hs94-stable20.inp
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ problem_id = hs94 # problem ID: basename of output filenames

<output0>
file_type = rst
dt = 1.0368E8
dt = 1.0368E6

<output1>
file_type = hst # History data dump
Expand All @@ -18,24 +18,24 @@ data_format = %15.12e
<output2>
file_type = netcdf # netcdf data dump
variable = prim # variables to be output
dt = 1.456E5 # time increment between outputs
dt = 3.456E5 # time increment between outputs

<output3>
file_type = netcdf # netcdf data dump
variable = uov # variables to be output
dt = 1.456E5 # time increment between outputs
dt = 3.456E5 # time increment between outputs

<output4>
file_type = dbg # netcdf data dump
dt = 1.456E2 # time increment between outputs
#<output4>
#file_type = dbg # netcdf data dump
#dt = 3.456E3 # time increment between outputs

<time>
cfl_number = 0.9 # The Courant, Friedrichs, & Lewy (CFL) Number
integrator = rk3 # integration method
xorder = 5 # horizontal reconstruction order
nlim = -1 # cycle limit
tlim = 2.0368E6 # time limit
shock_capture = true
tlim = 1000.0368E5 # time limit
#shock_capture = true

<mesh>
nx1 = 40 # Number of zones in X1-direction (radial)
Expand Down Expand Up @@ -67,7 +67,7 @@ nx3 = 20
<hydro>
gamma = 1.4 # gamma = C_p/C_v
grav_acc1 = -9.81 # gravity accelaration
implicit_flag = 1
implicit_flag = 9


<thermodynamics>
Expand All @@ -77,7 +77,7 @@ Rd = 287.
cp = 1004.

<problem>
Omega = 0. # 7.292E-5
Omega = 7.292E-5
Rp = 6.371E6
p0 = 1.E5 # surface pressure
Ts = 315. # surface temperature
Expand Down
24 changes: 23 additions & 1 deletion examples/2023-Chen-exo3/hs94.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -291,13 +291,16 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
}

void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
AllocateUserOutputVariables(6);
AllocateUserOutputVariables(9);
SetUserOutputVariableName(0, "temp");
SetUserOutputVariableName(1, "theta");
SetUserOutputVariableName(2, "lat");
SetUserOutputVariableName(3, "lon");
SetUserOutputVariableName(4, "vlat");
SetUserOutputVariableName(5, "vlon");
SetUserOutputVariableName(6, "Teq");
SetUserOutputVariableName(7, "Kv");
SetUserOutputVariableName(8, "Kt");
}

// \brif Output distributions of temperature and potential temperature.
Expand All @@ -321,5 +324,24 @@ void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
user_out_var(3, k, j, i) = lon;
user_out_var(4, k, j, i) = U;
user_out_var(5, k, j, i) = V;

Real kappa; // Rd/Cp
kappa = Rd / cp;

Real scaled_z = phydro->w(IPR, k, j, i) / p0;
Real Teq_p =
Ts - dT * _sqr(sin(lat)) - dtheta * log(scaled_z) * _sqr(cos(lat));
Teq_p *= pow(scaled_z, kappa);
Real Teq = (Teq_p > 200.) ? Teq_p : 200.;
user_out_var(6, k, j, i) = Teq;

Real sigma = phydro->w(IPR, k, j, i) / p0; // pmb->phydro->pbot(k,j);
Real sigma_p = (sigma - sigmab) / (1. - sigmab);
Real Kv = (sigma_p < 0.0) ? 0.0 : sigma_p * Kf;
user_out_var(7, k, j, i) = Kv;

sigma_p = (sigma_p < 0.0) ? 0.0 : sigma_p * _qur(cos(lat));
Real Kt = Ka + (Ks - Ka) * sigma_p;
user_out_var(8, k, j, i) = Kt;
}
}
5 changes: 5 additions & 0 deletions main.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,8 @@ int main(int argc, char **argv) {

std::uint64_t mbcnt = 0;
bool redo = false;
int max_redo = 1;
int attempts = 0;

while ((pmesh->time < pmesh->tlim) &&
(pmesh->nlim < 0 || pmesh->ncycle < pmesh->nlim)) {
Expand All @@ -139,6 +141,9 @@ int main(int argc, char **argv) {
pmesh->LoadAllStates();
pmesh->DecreaseTimeStep();
redo = true;
if (attempts++ > max_redo) {
throw RuntimeError("main", "Too many attempts to redo the step");
}
continue;
}

Expand Down
7 changes: 6 additions & 1 deletion src/snap/implicit/full_correction.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,13 @@ void ImplicitSolver::FullCorrection(AthenaArray<Real>& du,
// 0. forcing and volume matrix

Real gamma = pmy_block_->peos->GetGamma();
Real grav = pmy_block_->phydro->hsrc.GetG1();
Eigen::Matrix<Real, 5, 5> Phi, Dt, Bnds, Bnde, tmp;
Phi.setZero();
if (mydir_ == X1DIR) {
Phi(ivx, idn) = grav;
Phi(ien, ivx) = grav;
}

Dt.setIdentity();
Dt *= 1. / dt;
Expand Down Expand Up @@ -140,7 +146,6 @@ void ImplicitSolver::FullCorrection(AthenaArray<Real>& du,

// 5. fix boundary condition
if (first_block && !periodic_boundary) a[is] += b[is] * Bnds;

if (last_block && !periodic_boundary) a[ie] += c[ie] * Bnde;

// 6. solve tridiagonal system
Expand Down
55 changes: 36 additions & 19 deletions src/snap/implicit/solve_implicit_3d.cpp
Original file line number Diff line number Diff line change
@@ -1,15 +1,17 @@
// athena
#include <athena/mesh/mesh.hpp>
#include <athena/hydro/hydro.hpp>
#include <athena/stride_iterator.hpp>

// application
#include <application/exceptions.hpp>

// canoe
#include <configure.hpp>

// canoe
// exo3
#include <exo3/gnomonic_equiangle.hpp>
#include <exo3/cubed_sphere_utility.hpp>

// snap
#include "implicit_solver.hpp"
Expand All @@ -18,6 +20,10 @@
#include <glog/logging.h>
#endif

#ifdef CUBED_SPHERE
namespace cs = CubedSphereUtility;
#endif

void ImplicitSolver::SolveImplicit3D(AthenaArray<Real> &du, AthenaArray<Real> &w,
Real dt)
{
Expand Down Expand Up @@ -112,47 +118,58 @@ void ImplicitSolver::SolveImplicit3D(AthenaArray<Real> &du, AthenaArray<Real> &w
Real sin_theta = pco->GetSineCell(k, j);

for (int i = 0; i < w.GetDim1(); ++i) {
Real vy = w(IVY, k, j, i);
Real vz = w(IVZ, k, j, i);
w(IVY, k, j, i) += vz * cos_theta;
w(IVY, k, j, i) += w(IVZ, k, j, i) * cos_theta;
w(IVZ, k, j, i) *= sin_theta;

//cs::CovariantToContravariant(du_.at(k,j,i), cos_theta);
//du_(IVY, k, j, i) += du_(IVZ, k, j, i) * cos_theta;
//du_(IVZ, k, j, i) *= sin_theta;
}
#endif // CUBED_SPHERE

// do implicit
if (implicit_flag_ & (1 << 3))
if (implicit_flag_ & (1 << 3)) {
FullCorrection(du_, w, dt, k, j, is, ie);
else
} else {
PartialCorrection(du_, w, dt, k, j, is, ie);
}

#ifdef ENABLE_GLOG
// check for negative density and internal energy
for (int i = is; i <= ie; i++) {
LOG_IF(WARNING, ph->u(IEN,k,j,i) + du_(IEN,k,j,i) < 0.)
<< "rank = " << Globals::my_rank << ", (k,j,i) = "
<< "(" << k << "," << j << "," << i << ")"
<< ", u[IEN](before) = " << ph->u(IEN,k,j,i) + du(IEN,k,j,i)
<< ", u[IEN](after) = " << ph->u(IEN,k,j,i) + du_(IEN,k,j,i)
<< ", u[IVX](before) = " << ph->u(IVX,k,j,i) + du(IVX,k,j,i) << std::endl
<< ", u[IVX](after) = " << ph->u(IVX,k,j,i) + du_(IVX,k,j,i) << std::endl;
<< "(" << k << "," << j << "," << i << ")" << std::endl
<< "(before) u[IDN] = " << ph->u(IDN,k,j,i) + du(IDN,k,j,i)
<< ", u[IVX] = " << ph->u(IVX,k,j,i) + du(IVX,k,j,i)
<< ", u[IEN] = " << ph->u(IEN,k,j,i) + du(IEN,k,j,i) << std::endl

<< "(after) u[IDN] = " << ph->u(IDN,k,j,i) + du_(IDN,k,j,i)
<< ", u[IVX] = " << ph->u(IVX,k,j,i) + du_(IVX,k,j,i)
<< ", u[IEN] = " << ph->u(IEN,k,j,i) + du_(IEN,k,j,i) << std::endl;

LOG_IF(WARNING, ph->u(IDN,k,j,i) + du_(IDN,k,j,i) < 0.)
<< "rank = " << Globals::my_rank << ", (k,j,i) = "
<< "(" << k << "," << j << "," << i << ")"
<< ", u[IDN](before) = " << ph->u(IDN,k,j,i) + du(IDN,k,j,i)
<< ", u[IDN](after) = " << ph->u(IDN,k,j,i) + du_(IDN,k,j,i)
<< ", u[IVX](before) = " << ph->u(IVX,k,j,i) + du(IVX,k,j,i) << std::endl
<< ", u[IVX](after) = " << ph->u(IVX,k,j,i) + du_(IVX,k,j,i) << std::endl;
<< "(before) u[IDN] = " << ph->u(IDN,k,j,i) + du(IDN,k,j,i)
<< ", u[IVX] = " << ph->u(IVX,k,j,i) + du(IVX,k,j,i)
<< ", u[IEN] = " << ph->u(IEN,k,j,i) + du(IEN,k,j,i) << std::endl

<< "(after) u[IDN] = " << ph->u(IDN,k,j,i) + du_(IDN,k,j,i)
<< ", u[IVX] = " << ph->u(IVX,k,j,i) + du_(IVX,k,j,i)
<< ", u[IEN] = " << ph->u(IEN,k,j,i) + du_(IEN,k,j,i) << std::endl;
}
#endif // ENABLE_GLOG

// project back
// de-project velocity
#ifdef CUBED_SPHERE
for (int i = 0; i < w.GetDim1(); ++i) {
Real vy = w(IVY, k, j, i);
Real vz = w(IVZ, k, j, i);
w(IVY, k, j, i) -= vz / sin_theta * cos_theta;
w(IVZ, k, j, i) /= sin_theta;
w(IVY, k, j, i) -= w(IVZ, k, j, i) * cos_theta;

//du_(IVZ, k, j, i) /= sin_theta;
//du_(IVY, k, j, i) -= du_(IVZ, k, j) * cos_theta;
//cs::ContravariantToCovariant(du_.at(k, j, i), cos_theta);
}
#endif // CUBED_SPHERE
}
Expand Down
80 changes: 80 additions & 0 deletions tests/test_transform.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,80 @@
// C/C++
#include <algorithm>
#include <cmath>

// external
#include <gtest/gtest.h>

// athena
#include <athena/athena.hpp>

// exo3
#include <exo3/cubed_sphere_utility.hpp>

namespace cs = CubedSphereUtility;

TEST(CovariantToContravariant, test1) {
Real w[5] = {1., 2., 3., 4., 5.};
Real cth = 0.5;

cs::CovariantToContravariant(w, cth);

EXPECT_DOUBLE_EQ(w[0], 1.);
EXPECT_DOUBLE_EQ(w[1], 2.);
EXPECT_DOUBLE_EQ(w[2], 4./3.);
EXPECT_DOUBLE_EQ(w[3], 10./3.);
EXPECT_DOUBLE_EQ(w[4], 5.);
}

TEST(ContravariantToOrthogonal, test1) {
Real w[5] = {1., 2., 3., 4., 5.};
Real y[5] = {1., 2., 3., 4., 5.};
Real cth = 0.5;
Real sth = sqrt(1. - cth * cth);

cs::ContravariantToCovariant(y, cth);

Real ek1 = w[IVX] * y[IVX] + w[IVY] * y[IVY] + w[IVZ] * y[IVZ];

w[IVY] += w[IVZ] * cth;
w[IVZ] *= sth;

Real ek2 = w[IVX] * w[IVX] + w[IVY] * w[IVY] + w[IVZ] * w[IVZ];

EXPECT_DOUBLE_EQ(w[0], 1.);
EXPECT_DOUBLE_EQ(w[1], 2.);
EXPECT_DOUBLE_EQ(w[2], 5);
EXPECT_DOUBLE_EQ(w[3], sqrt(12));
EXPECT_DOUBLE_EQ(w[4], 5.);
EXPECT_DOUBLE_EQ(ek1, ek2);

w[IVZ] /= sth;
w[IVY] -= w[IVZ] * cth;

EXPECT_DOUBLE_EQ(w[0], 1.);
EXPECT_DOUBLE_EQ(w[1], 2.);
EXPECT_DOUBLE_EQ(w[2], 3.);
EXPECT_DOUBLE_EQ(w[3], 4.);
EXPECT_DOUBLE_EQ(w[4], 5.);
}

TEST(ContravariantToCovariant, test1) {
Real w[5] = {1., 2., 4./3., 10./3., 5.};
Real cth = 0.5;

cs::ContravariantToCovariant(w, cth);

EXPECT_DOUBLE_EQ(w[0], 1.);
EXPECT_DOUBLE_EQ(w[1], 2.);
EXPECT_DOUBLE_EQ(w[2], 3.);
EXPECT_DOUBLE_EQ(w[3], 4.);
EXPECT_DOUBLE_EQ(w[4], 5.);
}

int main(int argc, char **argv) {
testing::InitGoogleTest(&argc, argv);

int result = RUN_ALL_TESTS();

return result;
}

0 comments on commit 1d2a847

Please sign in to comment.