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

Cluster improvements and donut model #82

Merged
merged 27 commits into from
Sep 25, 2023
Merged

Cluster improvements and donut model #82

merged 27 commits into from
Sep 25, 2023

Conversation

pgrete
Copy link
Contributor

@pgrete pgrete commented Aug 18, 2023

Summary (should all be reflected in the updated doc)

  • Add multi-dimensional limited prolongation
  • Add tracer for jet material (from the kinetic jet launching region)
  • Add magic stellar heating
  • Allow to disable mass injection from the magnetic tower
  • Use persistent vector potential for tower (instead of alloc/free)
  • Allow to define char. lengthscale of initial perturbation (instead of normalized wavenumber)
  • Add donut model (and allow for runtime option to switch between)

Copy link
Contributor

@forrestglines forrestglines left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good -- I have some suggestions but none are blocking.

using parthenon::TE;
using parthenon::TopologicalElement;

struct ProlongateCellMinModMultiD {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we have an equation to reference for this? Stone 2020 Equation 5? With the caveat we're apply this to the conserveds?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, plus limited based on the AMReX implementation.
I added a pointer to both as a comment.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wrote that code, but I think it's correct nonetheless ;)

const int fj = (DIM > 1) ? (j - cjb.s) * 2 + jb.s : jb.s;
const int fk = (DIM > 2) ? (k - ckb.s) * 2 + kb.s : kb.s;

constexpr bool INCLUDE_X1 =
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So we include a dimension if:

  • If we're prolongating to a cell centered position
  • Dimensions along the plane of a face for a face element
  • Dimensions along the edge for an edge element

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, and to be honest, I adapted this logic from the Parthenon prolongation (without checking) trusting that it works there and we just needed limiting.

docs/cluster.md Outdated Show resolved Hide resolved
docs/cluster.md Outdated Show resolved Hide resolved
Comment on lines +72 to +83
// Normalize the thermal, kinetic, and magnetic mass fractions to sum to 1.0
if (enable_magnetic_tower_mass_injection_) {
thermal_mass_fraction_ = thermal_fraction_;
kinetic_mass_fraction_ = kinetic_fraction_;
magnetic_mass_fraction_ = magnetic_fraction_;
} else {
const auto total_mass_frac = thermal_fraction_ + kinetic_fraction_;
thermal_mass_fraction_ = thermal_fraction_ / total_mass_frac;
kinetic_mass_fraction_ = kinetic_fraction_ / total_mass_frac;
magnetic_mass_fraction_ = 0.0;
}

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In a commit I ultimately dropped, I had the mass fractions of each of the different mechanisms optionally decoupled from the energy fractions. Worth bringing back, would give a bit more functionality?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Potentially. The only reason I decided to go for this route was that it introduces only one additional parameter in the input file.
I suggest that we make further modifications (like explicit mass fraction parameters) when needed.

@@ -90,15 +86,16 @@ void MagneticTower::AddSrcTerm(parthenon::Real field_to_add, parthenon::Real mas
a_jb.e, a_ib.s, a_ib.e,
KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) {
// Compute and apply potential
auto &A = A_pack(b);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Any particular reason to make this switch? Maybe caching that first array access into the block?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Which switch specifically?
Previously, the potential was an explicit ParArray5D allocated an deallocated on every single call.
Now it's a standard Parthenon field that we access through the packing machinery.

Or are you referring to not using the b index and access A_pack(b,n,k,j,i) in the following lines?
In that case, I made the switch because we typically extract the block pointer also for prim and cons fields.

a_r = 0.0;
a_theta = field_ * l_scale_ * (r / l_scale_) * exp_r2_h2;
a_h = field_ * l_scale_ * alpha_ / 2.0 * exp_r2_h2;
if (potential_ == MagneticTowerPotential::donut) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure I like this if called everytime, I'd like to optimize this out at compiler time and an if statement elsewhere.

Not necessary for this PR, however

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm am sure that I don't like it ;)
I just went for it for simplicity now.
I started to write a template MagneticTowerObj but then realized that we call the constructor at three places and I didn't want to add the identical logic at three places (or introduce std::variant or similar).
Therefore, I decided to keep it simple and revisit once we add even more options.

src/pgen/cluster/stellar_feedback.cpp Show resolved Hide resolved
src/pgen/cluster/stellar_feedback.cpp Outdated Show resolved Hide resolved
@forrestglines
Copy link
Contributor

Merged in the cluster clips/stellar feedback reductions and cold gas/AGN extent reductions into this PR

@pgrete
Copy link
Contributor Author

pgrete commented Sep 20, 2023

I added a couple of minor fixes (param names, prim scalar tracer for the extent (rather than the scalar density, and formatting, see 2bdf1f8...36d5139

I think this is ready for a merge (as the PR is growing, I suggest to make further modifications separately).
Any objections @forrestglines ?

@forrestglines
Copy link
Contributor

I put up one last commit that just removes commented out code for the alternative methods to include cooling.

I wanted to test this on Frontier one last time from a clean build to make sure everything is working. I got distracted trying to get ROCM 5.6.0 and ROCM 5.7.0 working but the ROCM 5.4.3 build works as expected. Merging this in.

@forrestglines forrestglines merged commit 6dc42ec into main Sep 25, 2023
4 checks passed
@pgrete pgrete deleted the pgrete/donut branch March 18, 2024 10:44
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants