Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fully-implicit cosmic-ray diffusion module using Multigrid #550

Merged
merged 35 commits into from
Mar 12, 2024
Merged
Changes from 1 commit
Commits
Show all changes
35 commits
Select commit Hold shift + click to select a range
2cfce54
Added a flag to control red-black iteration or not.
tomidakn Nov 19, 2023
bf23edc
**incomplete commit** implementing MG CR Diffusion
tomidakn Nov 20, 2023
7330160
*** incomplete commit *** working on multigrid CR diffusion
tomidakn Nov 21, 2023
a038ba0
*** incomplete commit *** implemented RestrictCoefficients
tomidakn Nov 22, 2023
e895429
*** incomplete commit ***
tomidakn Nov 23, 2023
3e85e54
*** incomplete commit *** changed MG coefficients to AthenaArray.
tomidakn Nov 23, 2023
cea80ab
*** incomplete commit *** save
tomidakn Nov 24, 2023
5a7b26c
*** incomplete commit *** save
tomidakn Nov 26, 2023
0f958a6
*** incomplete commit *** Fixed a bug related to ExpandPhysicalBounda…
tomidakn Nov 26, 2023
50eda93
*** incomplete commit *** implemented matrix element calculation
tomidakn Nov 26, 2023
a19d032
First try - it does not converge.
tomidakn Nov 26, 2023
4468e1f
*** incomplete commit *** fixed a bug in FASRHS
tomidakn Nov 27, 2023
dd20fce
Seems working but need to improve convergence.
tomidakn Nov 27, 2023
b82fc5d
Simplified OpenMP (not tested yet)
tomidakn Nov 28, 2023
b5a271c
Implemented Jacobi-RB smoother
tomidakn Nov 30, 2023
8425eb7
Seems to be working! Still need to test AMR.
tomidakn Nov 30, 2023
33de17a
Implemented coefficient boundary and mask functions
tomidakn Dec 2, 2023
339c8dc
*** incomplete commit *** First implemention of boundary communicatio…
tomidakn Dec 2, 2023
f22fdd8
MPI works!
tomidakn Dec 3, 2023
ab4cba1
Fixed a bug in memory allocation
tomidakn Dec 5, 2023
258a54d
Fixed a bug in MG mask functions
tomidakn Dec 5, 2023
291ffc7
Fixed a memory allocation bug in AMR.
tomidakn Dec 7, 2023
6e07efb
Added some new features on Multigrid
tomidakn Dec 13, 2023
3e59ebc
Fixing style errors
tomidakn Dec 24, 2023
3ae2638
Fixed inconsistent class/structure definitions
tomidakn Dec 24, 2023
148a6e6
Fix unused variables
tomidakn Dec 25, 2023
6e984c6
Added CR diffusion in the pgen_compile.py
tomidakn Dec 25, 2023
34f75f0
Added a line in pgen_compile.py
tomidakn Dec 25, 2023
c747d3c
Fixed a bug in pgen_compile.py
tomidakn Dec 25, 2023
61734ba
Merge branch 'master' into crdiff
felker Dec 28, 2023
b19bfe1
Fixed a merging error
tomidakn Dec 30, 2023
24810d5
Fixed an error in pgen_compile.py
tomidakn Dec 30, 2023
2268de9
Removed zeta from CR diffusion, and minor fixes
tomidakn Mar 2, 2024
dbfad59
Merge branch 'master' of https://github.com/PrincetonUniversity/athen…
tomidakn Mar 2, 2024
640a7b4
omega for self-gravity can be set in the input file
tomidakn Mar 2, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Seems to be working! Still need to test AMR.
tomidakn committed Nov 30, 2023
commit 8425eb719b706223b7dd8f93ba9d551ceb26d518
43 changes: 24 additions & 19 deletions inputs/cosmic_ray/athinput.cr_diffusion_mg
Original file line number Diff line number Diff line change
@@ -13,11 +13,11 @@ dt = 0.01 # time increment between outputs
<output1>
file_type = vtk # Binary data dump
variable = prim # variables to be output
dt = 0.01 # time increment between outputs
dcycle = 1 # time increment between outputs

<time>
cfl_number = 0.3 # The Courant, Friedrichs, & Lewy (CFL) Number
nlim = 0 # cycle limit
nlim = 1 # cycle limit
tlim = 1.0 # time limit
integrator = vl2 # time integration algorithm
xorder = 2 # order of spatial reconstruction
@@ -27,20 +27,20 @@ ncycle_out = 1 # interval for stdout summary info
nx1 = 64 # Number of zones in X1-direction
x1min = -0.5 # minimum value of X1
x1max = 0.5 # maximum value of X1
ix1_bc = periodic # inner-X1 boundary flag
ox1_bc = periodic # outer-X1 boundary flag
ix1_bc = outflow # inner-X1 boundary flag
ox1_bc = outflow # outer-X1 boundary flag

nx2 = 64 # Number of zones in X2-direction
x2min = -0.5 # minimum value of X2
x2max = 0.5 # maximum value of X2
ix2_bc = periodic # inner-X2 boundary flag
ox2_bc = periodic # outer-X2 boundary flag
ix2_bc = outflow # inner-X2 boundary flag
ox2_bc = outflow # outer-X2 boundary flag

nx3 = 64 # Number of zones in X3-direction
x3min = -0.5 # minimum value of X3
x3max = 0.5 # maximum value of X3
ix3_bc = periodic # inner-X3 boundary flag
ox3_bc = periodic # outer-X3 boundary flag
ix3_bc = outflow # inner-X3 boundary flag
ox3_bc = outflow # outer-X3 boundary flag

<meshblock>
nx1 = 64
@@ -53,18 +53,23 @@ iso_sound_speed = 0.4082482905 # equivalent to sqrt(gamma*p/d) for p=0.1, d=1

<crdiffusion>
mgmode = FMG
niteration = 1
#threshold = 0.0
fas = true
npresmooth = 2
npostsmooth = 2
omega = 1.3
niteration = 10
#threshold = 0.00001
output_defect = true
ix1_bc = periodic
ox1_bc = periodic
ix2_bc = periodic
ox2_bc = periodic
ix3_bc = periodic
ox3_bc = periodic
ix1_bc = user
ox1_bc = user
ix2_bc = user
ox2_bc = user
ix3_bc = zerograd
ox3_bc = zerograd

Dpara = 1.0
Dperp = 0.01
Lambda = 0.0
Dpara = 100.0
Dperp = 1.0
Lambda = 100.0
zeta_factor = 1.0

<problem>
7 changes: 3 additions & 4 deletions src/crdiffusion/crdiffusion.cpp
Original file line number Diff line number Diff line change
@@ -45,6 +45,7 @@ CRDiffusion::CRDiffusion(MeshBlock *pmb, ParameterInput *pin) :
Dpara_ = pin->GetReal("crdiffusion", "Dpara");
Dperp_ = pin->GetReal("crdiffusion", "Dperp");
Lambda_ = pin->GetReal("crdiffusion", "Lambda");
zeta_factor_ = pin->GetReal("crdiffusion", "zeta_factor");

output_defect = pin->GetOrAddBoolean("crdiffusion", "output_defect", false);
if (output_defect)
@@ -84,8 +85,7 @@ void CRDiffusion::CalculateCoefficients(const AthenaArray<Real> &w,
jl -= NGHOST, ju += NGHOST;
if (pmy_block->pmy_mesh->f3)
kl -= NGHOST, ku += NGHOST;
constexpr Real Dunit = 1.0, nlunit = 1.0;
Real Dpara = Dpara_ * Dunit, Dperp = Dperp_ * Dunit, Lambda = Lambda_ * nlunit;
Real Dpara = Dpara_, Dperp = Dperp_, Lambda = Lambda_;

if (MAGNETIC_FIELDS_ENABLED) {
for (int k = kl; k <= ku; ++k) {
@@ -126,7 +126,6 @@ void CRDiffusion::CalculateCoefficients(const AthenaArray<Real> &w,
//! \fn void CRDiffusion::CalculateIonizationRate(const AthenaArray<Real> &w)
//! \brief Calculate Ionization Rate from the Cosmic-ray density
void CRDiffusion::CalculateIonizationRate(const AthenaArray<Real> &w) {
Real zeta_factor = 1.0;
int il = pmy_block->is - NGHOST, iu = pmy_block->ie + NGHOST;
int jl = pmy_block->js, ju = pmy_block->je;
int kl = pmy_block->ks, ku = pmy_block->ke;
@@ -137,7 +136,7 @@ void CRDiffusion::CalculateIonizationRate(const AthenaArray<Real> &w) {
for (int k = kl; k <= ku; ++k) {
for (int j = jl; j <= ju; ++j) {
for (int i = il; i <= iu; ++i)
zeta(k, j, i) = zeta_factor * w(IDN, k, j, i) * ecr(k, j, i);
zeta(k, j, i) = zeta_factor_ * w(IDN, k, j, i) * ecr(k, j, i);
}
}

2 changes: 1 addition & 1 deletion src/crdiffusion/crdiffusion.hpp
Original file line number Diff line number Diff line change
@@ -57,7 +57,7 @@ class CRDiffusion {

private:
int refinement_idx_;
Real Dpara_, Dperp_, Lambda_;
Real Dpara_, Dperp_, Lambda_, zeta_factor_;
};

#endif // CRDIFFUSION_CRDIFFUSION_HPP_
145 changes: 112 additions & 33 deletions src/crdiffusion/mg_crdiffusion.cpp
Original file line number Diff line number Diff line change
@@ -50,7 +50,23 @@ MGCRDiffusionDriver::MGCRDiffusionDriver(Mesh *pm, ParameterInput *pin)
niter_ = pin->GetOrAddInteger("crdiffusion", "niteration", -1);
ffas_ = pin->GetOrAddBoolean("crdiffusion", "fas", ffas_);
omega_ = pin->GetOrAddReal("crdiffusion", "omega", omega_);
redblack_ = true;
npresmooth_ = pin->GetOrAddReal("crdiffusion", "npresmooth", npresmooth_);
npostsmooth_ = pin->GetOrAddReal("crdiffusion", "npostsmooth", npostsmooth_);
std::string smoother = pin->GetOrAddString("crdiffusion", "smoother", "jacobi-rb");
if (smoother == "jacobi-rb") {
fsmoother_ = 1;
redblack_ = true;
} else if (smoother == "jacobi-double") {
fsmoother_ = 0;
redblack_ = true;
}else { // jacobi
fsmoother_ = 0;
redblack_ = false;
}
std::string prol = pin->GetOrAddString("crdiffusion", "prolongation", "trilinear");
if (prol == "tricubic")
fprolongation_ = 1;

std::string m = pin->GetOrAddString("crdiffusion", "mgmode", "none");
std::transform(m.begin(), m.end(), m.begin(), ::tolower);
if (m == "fmg") {
@@ -125,7 +141,7 @@ MGCRDiffusionDriver::~MGCRDiffusionDriver() {
//! \brief MGCRDiffusion constructor

MGCRDiffusion::MGCRDiffusion(MGCRDiffusionDriver *pmd, MeshBlock *pmb)
: Multigrid(pmd, pmb, 1), omega_(pmd->omega_) {
: Multigrid(pmd, pmb, 1), omega_(pmd->omega_), fsmoother_(pmd->fsmoother_) {
btype = btypef = BoundaryQuantity::mg;
pmgbval = new MGBoundaryValues(this, mg_block_bcs_);
}
@@ -145,7 +161,6 @@ MGCRDiffusion::~MGCRDiffusion() {
//! \brief load the data and solve

void MGCRDiffusionDriver::Solve(int stage, Real dt) {
std::cout << "Solve" << std::endl;
// Construct the Multigrid array
vmg_.clear();
for (int i = 0; i < pmy_mesh_->nblocal; ++i)
@@ -216,12 +231,46 @@ void MGCRDiffusion::Smooth(AthenaArray<Real> &u, const AthenaArray<Real> &src,
Real dx2 = SQR(dx);
Real isix = omega_/6.0;
color ^= pmy_driver_->coffset_;

if (th == true && (ku-kl) >= minth_) {
AthenaArray<Real> &work = temp[0];
if (fsmoother_ == 1) { // jacobi-rb
if (th == true && (ku-kl) >= minth_) {
AthenaArray<Real> &work = temp[0];
#pragma omp parallel num_threads(pmy_driver_->nthreads_)
{
{
#pragma omp for
for (int k=kl; k<=ku; k++) {
for (int j=jl; j<=ju; j++) {
int c = (color + k + j) & 1;
#pragma ivdep
for (int i=il+c; i<=iu; i+=2) {
Real M = matrix(CCM,k,j,i)*u(k,j,i-1) + matrix(CCP,k,j,i)*u(k,j,i+1)
+ matrix(CMC,k,j,i)*u(k,j-1,i) + matrix(CPC,k,j,i)*u(k,j+1,i)
+ matrix(MCC,k,j,i)*u(k-1,j,i) + matrix(PCC,k,j,i)*u(k+1,j,i)
+ matrix(CMM,k,j,i)*u(k,j-1,i-1) + matrix(CMP,k,j,i)*u(k,j-1,i+1)
+ matrix(CPM,k,j,i)*u(k,j+1,i-1) + matrix(CPP,k,j,i)*u(k,j+1,i+1)
+ matrix(MCM,k,j,i)*u(k-1,j,i-1) + matrix(MCP,k,j,i)*u(k-1,j,i+1)
+ matrix(PCM,k,j,i)*u(k+1,j,i-1) + matrix(PCP,k,j,i)*u(k+1,j,i+1)
+ matrix(MMC,k,j,i)*u(k-1,j-1,i) + matrix(MPC,k,j,i)*u(k-1,j+1,i)
+ matrix(PMC,k,j,i)*u(k+1,j-1,i) + matrix(PPC,k,j,i)*u(k+1,j+1,i);
work(k,j,i) = (src(k,j,i) - M) / matrix(CCC,k,j,i);
}
}
}
#pragma omp for
for (int k=kl; k<=ku; k++) {
for (int j=jl; j<=ju; j++) {
int c = (color + k + j) & 1;
#pragma ivdep
for (int i=il+c; i<=iu; i+=2)
u(k,j,i) += omega_ * (work(k,j,i) - u(k,j,i));
}
}
}
} else {
int t = 0;
#ifdef OPENMP_PARALLEL
t = omp_get_thread_num();
#endif
AthenaArray<Real> &work = temp[t];
for (int k=kl; k<=ku; k++) {
for (int j=jl; j<=ju; j++) {
int c = (color + k + j) & 1;
@@ -236,11 +285,10 @@ void MGCRDiffusion::Smooth(AthenaArray<Real> &u, const AthenaArray<Real> &src,
+ matrix(PCM,k,j,i)*u(k+1,j,i-1) + matrix(PCP,k,j,i)*u(k+1,j,i+1)
+ matrix(MMC,k,j,i)*u(k-1,j-1,i) + matrix(MPC,k,j,i)*u(k-1,j+1,i)
+ matrix(PMC,k,j,i)*u(k+1,j-1,i) + matrix(PPC,k,j,i)*u(k+1,j+1,i);
work(k,j,i) = (src(k,j,i) - M) / matrix(CCC,k,j,i);
work(k,j,i) = (src(k,j,i) - M) / matrix(CCC,k,j,i);
}
}
}
#pragma omp for
for (int k=kl; k<=ku; k++) {
for (int j=jl; j<=ju; j++) {
int c = (color + k + j) & 1;
@@ -250,36 +298,67 @@ void MGCRDiffusion::Smooth(AthenaArray<Real> &u, const AthenaArray<Real> &src,
}
}
}
} else {
int t = 0;
} else { // jacobi
if (th == true && (ku-kl) >= minth_) {
AthenaArray<Real> &work = temp[0];
#pragma omp parallel num_threads(pmy_driver_->nthreads_)
{
#pragma omp for
for (int k=kl; k<=ku; k++) {
for (int j=jl; j<=ju; j++) {
#pragma ivdep
for (int i=il; i<=iu; i++) {
Real M = matrix(CCM,k,j,i)*u(k,j,i-1) + matrix(CCP,k,j,i)*u(k,j,i+1)
+ matrix(CMC,k,j,i)*u(k,j-1,i) + matrix(CPC,k,j,i)*u(k,j+1,i)
+ matrix(MCC,k,j,i)*u(k-1,j,i) + matrix(PCC,k,j,i)*u(k+1,j,i)
+ matrix(CMM,k,j,i)*u(k,j-1,i-1) + matrix(CMP,k,j,i)*u(k,j-1,i+1)
+ matrix(CPM,k,j,i)*u(k,j+1,i-1) + matrix(CPP,k,j,i)*u(k,j+1,i+1)
+ matrix(MCM,k,j,i)*u(k-1,j,i-1) + matrix(MCP,k,j,i)*u(k-1,j,i+1)
+ matrix(PCM,k,j,i)*u(k+1,j,i-1) + matrix(PCP,k,j,i)*u(k+1,j,i+1)
+ matrix(MMC,k,j,i)*u(k-1,j-1,i) + matrix(MPC,k,j,i)*u(k-1,j+1,i)
+ matrix(PMC,k,j,i)*u(k+1,j-1,i) + matrix(PPC,k,j,i)*u(k+1,j+1,i);
work(k,j,i) = (src(k,j,i) - M) / matrix(CCC,k,j,i);
}
}
}
#pragma omp for
for (int k=kl; k<=ku; k++) {
for (int j=jl; j<=ju; j++) {
#pragma ivdep
for (int i=il; i<=iu; i++)
u(k,j,i) += omega_ * (work(k,j,i) - u(k,j,i));
}
}
}
} else {
int t = 0;
#ifdef OPENMP_PARALLEL
t = omp_get_thread_num();
t = omp_get_thread_num();
#endif
AthenaArray<Real> &work = temp[t];
for (int k=kl; k<=ku; k++) {
for (int j=jl; j<=ju; j++) {
int c = (color + k + j) & 1;
AthenaArray<Real> &work = temp[t];
for (int k=kl; k<=ku; k++) {
for (int j=jl; j<=ju; j++) {
#pragma ivdep
for (int i=il+c; i<=iu; i+=2) {
Real M = matrix(CCM,k,j,i)*u(k,j,i-1) + matrix(CCP,k,j,i)*u(k,j,i+1)
+ matrix(CMC,k,j,i)*u(k,j-1,i) + matrix(CPC,k,j,i)*u(k,j+1,i)
+ matrix(MCC,k,j,i)*u(k-1,j,i) + matrix(PCC,k,j,i)*u(k+1,j,i)
+ matrix(CMM,k,j,i)*u(k,j-1,i-1) + matrix(CMP,k,j,i)*u(k,j-1,i+1)
+ matrix(CPM,k,j,i)*u(k,j+1,i-1) + matrix(CPP,k,j,i)*u(k,j+1,i+1)
+ matrix(MCM,k,j,i)*u(k-1,j,i-1) + matrix(MCP,k,j,i)*u(k-1,j,i+1)
+ matrix(PCM,k,j,i)*u(k+1,j,i-1) + matrix(PCP,k,j,i)*u(k+1,j,i+1)
+ matrix(MMC,k,j,i)*u(k-1,j-1,i) + matrix(MPC,k,j,i)*u(k-1,j+1,i)
+ matrix(PMC,k,j,i)*u(k+1,j-1,i) + matrix(PPC,k,j,i)*u(k+1,j+1,i);
work(k,j,i) = (src(k,j,i) - M) / matrix(CCC,k,j,i);
for (int i=il; i<=iu; i++) {
Real M = matrix(CCM,k,j,i)*u(k,j,i-1) + matrix(CCP,k,j,i)*u(k,j,i+1)
+ matrix(CMC,k,j,i)*u(k,j-1,i) + matrix(CPC,k,j,i)*u(k,j+1,i)
+ matrix(MCC,k,j,i)*u(k-1,j,i) + matrix(PCC,k,j,i)*u(k+1,j,i)
+ matrix(CMM,k,j,i)*u(k,j-1,i-1) + matrix(CMP,k,j,i)*u(k,j-1,i+1)
+ matrix(CPM,k,j,i)*u(k,j+1,i-1) + matrix(CPP,k,j,i)*u(k,j+1,i+1)
+ matrix(MCM,k,j,i)*u(k-1,j,i-1) + matrix(MCP,k,j,i)*u(k-1,j,i+1)
+ matrix(PCM,k,j,i)*u(k+1,j,i-1) + matrix(PCP,k,j,i)*u(k+1,j,i+1)
+ matrix(MMC,k,j,i)*u(k-1,j-1,i) + matrix(MPC,k,j,i)*u(k-1,j+1,i)
+ matrix(PMC,k,j,i)*u(k+1,j-1,i) + matrix(PPC,k,j,i)*u(k+1,j+1,i);
work(k,j,i) = (src(k,j,i) - M) / matrix(CCC,k,j,i);
}
}
}
}
for (int k=kl; k<=ku; k++) {
for (int j=jl; j<=ju; j++) {
int c = (color + k + j) & 1;
for (int k=kl; k<=ku; k++) {
for (int j=jl; j<=ju; j++) {
#pragma ivdep
for (int i=il+c; i<=iu; i+=2)
u(k,j,i) += omega_ * (work(k,j,i) - u(k,j,i));
for (int i=il; i<=iu; i++)
u(k,j,i) += omega_ * (work(k,j,i) - u(k,j,i));
}
}
}
}
2 changes: 2 additions & 0 deletions src/crdiffusion/mg_crdiffusion.hpp
Original file line number Diff line number Diff line change
@@ -50,6 +50,7 @@ class MGCRDiffusion : public Multigrid {

private:
Real omega_;
int fsmoother_;
};


@@ -67,6 +68,7 @@ class MGCRDiffusionDriver : public MultigridDriver {
private:
CRDiffusionBoundaryTaskList *crtlist_;
Real omega_;
int fsmoother_;
};

#endif // CRDIFFUSION_MG_CRDIFFUSION_HPP_
2 changes: 2 additions & 0 deletions src/gravity/mg_gravity.cpp
Original file line number Diff line number Diff line change
@@ -45,6 +45,8 @@ MGGravityDriver::MGGravityDriver(Mesh *pm, ParameterInput *pin)
eps_ = pin->GetOrAddReal("gravity", "threshold", -1.0);
niter_ = pin->GetOrAddInteger("gravity", "niteration", -1);
ffas_ = pin->GetOrAddBoolean("gravity", "fas", ffas_);
npresmooth_ = pin->GetOrAddReal("gravity", "npresmooth", npresmooth_);
npostsmooth_ = pin->GetOrAddReal("gravity", "npostsmooth", npostsmooth_);
redblack_ = true;
std::string m = pin->GetOrAddString("gravity", "mgmode", "none");
std::transform(m.begin(), m.end(), m.begin(), ::tolower);
Loading