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

New general anelastic non-dimensional subroutine (reference_type = 5) #450

Merged
merged 14 commits into from
Jun 14, 2023

Conversation

illorenzo7
Copy link
Contributor

@illorenzo7 illorenzo7 commented Mar 17, 2023

Implements to subroutine

Polytropic_Reference_ND_General()

in PDE_Coefficients. Note that my other pull request (#447) is included here. This new functionality allows user to specify simulation using Prandtl number, Rayleigh number (flux or otherwise), Ekman number, magnetic Prandtl number, and four number characterizing the polytrop ice: specific-heat ratio (default 5/3), poly tropic index, number of density scale heights, and shell aspect ratio.

User can choose between non-dimensionalizing the reference state using the profile values at the inner radius, outer radius, or volume-averaged over the shell (default).

User can non-dimensionalize time using either the rotation rate (time = 1/Omega_0) or viscous diffusion time.

User can non-dimensionalize the length scale using an arbitrary value, with default being the shell depth.

Endif

! Determine how user wants to specify Ra, Chi_A_Ra, and B_visc.
! (The modified versions or not).
Copy link
Contributor

Choose a reason for hiding this comment

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

We should have a branch for a model without rotation.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Actually I think now it's fine. At most, we could have the code throw a warning if they set rotation = .false but decide to set a non-zero Modified_Rayleigh_Number (say). But we don't do anything like that for e.g. the Magnetic_Prandtl_Number (the user can set that with magnetism=.false. and no warning will be thrown).


denom1 = (1.0d0 - exp(-poly_Nrho/poly_n))/(1.d0 - Aspect_Ratio)
denom2 = -(Aspect_Ratio - exp(-poly_Nrho/poly_n))/(Aspect_Ratio*Shell_Depth)
nsquared = (rmin/radius)**3 / (denom1 + denom2*radius)
Copy link
Contributor

Choose a reason for hiding this comment

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

Could we explicily write dsdr in terms of dlnrho and dlnP or dlnT and then multiply by the appropriate amplitude control parameter?

Call Integrate_in_radius(nsquared,norm)
norm = four_pi*norm/shell_volume
nsquared = nsquared/norm

Copy link
Contributor

Choose a reason for hiding this comment

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

Can you compute dissipation number in a similar way to the new way of calculating dsdr?

Copy link
Contributor Author

@illorenzo7 illorenzo7 Jun 13, 2023

Choose a reason for hiding this comment

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

I think what I just did is the cleanest way to do it (no need to integrate temperature and gravity a second time).

Copy link
Contributor

@feathern feathern 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. Thanks for going over this with me!

@feathern feathern merged commit 70adb26 into geodynamics:main Jun 14, 2023
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.

2 participants