Skip to content

Add curl momentum diagnostics #625

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.

Sign up for GitHub

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

Draft
Ayesha714 wants to merge 15 commits into
base: main
Choose a base branch
from
Draft

Add curl momentum diagnostics #625

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
09bc5e3
Add curl momentum diagnostics
Jul 17, 2026
19c9e4e
Removing line
cianwilson Jul 17, 2026
7db7fed
Adding a placeholder file to contain the curl of the momentum equatio…
cianwilson Jul 17, 2026
6c85b5f
Update curl momentum equation diagnostics
Jul 17, 2026
a926426
Refactor file (untested).
cianwilson Jul 17, 2026
1af6c41
Merge branch 'main' into Curl_Momentum_Diagnostics
cianwilson Jul 17, 2026
4276c6a
Merge branch geodynamics:main into Curl_Momentum_Diagnostics
cianwilson Jul 17, 2026
e974fbe
Implementation of the derivatives of the viscous forces. Untested bu…
cianwilson Jul 19, 2026
ceaa224
Trying to update docs to reflect new curl of momentum quantity codes.…
cianwilson Jul 19, 2026
192d1ca
Adding an untested example of turning the calculated gradients of the…
cianwilson Jul 19, 2026
ce21898
Remove comments that have been resolved by #633
cianwilson Jul 22, 2026
3d31fa0
Question about using new curl momentum diagnostics module
Jul 22, 2026
5382d47
Fix for the failing build test.
cianwilson Jul 22, 2026
07b8cd3
A logic fix for the squared versions of the curl of the viscous forces.
cianwilson Jul 22, 2026
aa96b56
A series of fixes for the curl of momentum diagnostics.
cianwilson Jul 24, 2026
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
2 changes: 2 additions & 0 deletions doc/build_menu_pdf.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,7 @@ def parse_offset_line(codeline):
inprefs.append(['current_density' ,2])
inprefs.append(['magnetic_energy' ,2])
inprefs.append(['momentum_equation' ,2])
inprefs.append(['curl_momentum_equation' ,2])
inprefs.append(['thermal_equation' ,2])

inprefs.append(['induction_equation' ,2])
Expand All @@ -130,6 +131,7 @@ def parse_offset_line(codeline):
page_titles.append('Current Density')
page_titles.append('Magnetic Energy')
page_titles.append('Momentum Equation')
page_titles.append('Curl Momentum Equation')
page_titles.append('Thermal Energy Equation')
page_titles.append('Induction Equation')
page_titles.append('Angular Momentum Equation')
Expand Down
85 changes: 85 additions & 0 deletions doc/source/diagnostic_codes/curl_momentum_equation.rst
View file
Open in desktop

Large diffs are not rendered by default.

1 change: 1 addition & 0 deletions doc/source/diagnostic_codes/qcodes.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ Output Quantity Codes
current_density
magnetic_energy
momentum_equation
curl_momentum_equation
thermal_equation
induction_equation
amom_equation
Expand Down
16 changes: 15 additions & 1 deletion src/Diagnostics/Diagnostics_Base.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,7 @@ Module Diagnostics_Base
include "magnetic_energy_codes.F"

include "momentum_equation_codes.F"
include "curl_momentum_equation_codes.F"
include "thermal_equation_codes.F"
include "induction_equation_codes.F"

Expand Down Expand Up @@ -131,7 +132,6 @@ Module Diagnostics_Base

Logical :: need_second_derivatives = .false.


!////////////////////////////////////////////////////////////////////////////
! Variables related to mean-correction
! (we only correct radial terms, but retain logic for horizontal terms)
Expand All @@ -152,6 +152,20 @@ Module Diagnostics_Base
Integer :: vfp_r, vfp_t, vfp_p
Integer :: vfm_r, vfm_t, vfm_p

! A special buffer used for holding first derivatives of the viscous forces at output time
Type(SphericalBuffer) :: d_vforce_buffer
Integer :: dvf_r_dt, dvf_r_dp
Integer :: dvf_t_dr, dvf_t_dp
Integer :: dvf_p_dr, dvf_p_dt
Integer :: dvfp_r_dt, dvfp_r_dp
Integer :: dvfp_t_dr, dvfp_t_dp
Integer :: dvfp_p_dr, dvfp_p_dt
Integer :: dvfm_r_dt, dvfm_r_dp
Integer :: dvfm_t_dr, dvfm_t_dp
Integer :: dvfm_p_dr, dvfm_p_dt

Logical :: need_vforce_derivatives = .false.

Contains


Expand Down
805 changes: 805 additions & 0 deletions src/Diagnostics/Diagnostics_Curl_Momentum.F90
View file
Open in desktop

Large diffs are not rendered by default.

12 changes: 10 additions & 2 deletions src/Diagnostics/Diagnostics_Interface.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@ Module Diagnostics_Interface
Use Diagnostics_Base

Use Diagnostics_Second_Derivatives
Use Diagnostics_Curl_Momentum

Use Diagnostics_Mean_Correction

Expand Down Expand Up @@ -174,7 +175,6 @@ Subroutine PS_Output(buffer,iteration, current_time)
Call ComputeM0(buffer,m0_values)
Call Compute_Fluctuations(buffer)

Call Initialize_Viscous_Force()
Call Initialize_Mean_Correction()


Expand All @@ -189,6 +189,10 @@ Subroutine PS_Output(buffer,iteration, current_time)
over_n_phi = 1.0d0/dble(n_phi)

Call Viscous_Force(buffer) ! Pre-calculate the viscous forces and place them in the vforce_buffer

if (need_vforce_derivatives) then
Call Grad_Viscous_Force()
endif

Call Mean_Correction(buffer) ! Remove ell=0 component from radial and theta forces

Expand Down Expand Up @@ -220,6 +224,7 @@ Subroutine PS_Output(buffer,iteration, current_time)
Call Compute_Angular_Momentum_Balance(buffer)
Call Compute_Inertial_Terms(buffer)
Call Compute_Linear_Forces(buffer)
Call Compute_Curl_Momentum_Forces(buffer)

Call Compute_KE_Flux(buffer)

Expand Down Expand Up @@ -257,7 +262,6 @@ Subroutine PS_Output(buffer,iteration, current_time)
Call d2buffer%deconstruct('p3a')
DeAllocate(d2_ell0,d2_m0,d2_fbuffer)
ENDIF
Call Finalize_Viscous_Force()
Call Finalize_Mean_Correction()
Endif ! time_to_output(iteration)
End Subroutine PS_Output
Expand Down Expand Up @@ -338,6 +342,10 @@ Subroutine Initialize_Diagnostics()

Call Initialize_Second_Derivatives()

Call Initialize_Viscous_Force()

Call Initialize_Grad_Viscous_Force()

Call Initialize_Diagnostics_Buffer()

End Subroutine Initialize_Diagnostics
Expand Down
111 changes: 76 additions & 35 deletions src/Diagnostics/Diagnostics_Linear_Forces.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -517,58 +517,82 @@ Subroutine Initialize_Viscous_Force()
Implicit None
nvf = 0

If (compute_quantity(viscous_force_r) .or. &
compute_quantity(visc_work) .or. &
compute_quantity(viscous_mforce_r)) Then
If (sometimes_compute(viscous_force_r) .or. &
sometimes_compute(visc_work) .or. &
sometimes_compute(curl_viscous_force_theta) .or. &
sometimes_compute(curl_viscous_force_theta_squared) .or. &
sometimes_compute(curl_viscous_force_phi) .or. &
sometimes_compute(curl_viscous_force_phi_squared) .or. &
sometimes_compute(viscous_mforce_r)) Then
nvf = nvf+1
vf_r = nvf
Endif

If (compute_quantity(viscous_force_theta) .or. &
compute_quantity(visc_work)) Then
If (sometimes_compute(viscous_force_theta) .or. &
sometimes_compute(visc_work) .or. &
sometimes_compute(curl_viscous_force_r) .or. &
sometimes_compute(curl_viscous_force_r_squared) .or. &
sometimes_compute(curl_viscous_force_phi) .or. &
sometimes_compute(curl_viscous_force_phi_squared)) Then
nvf = nvf + 1
vf_t = nvf
Endif

If (compute_quantity(viscous_force_phi) .or. &
compute_quantity(visc_work)) Then
If (sometimes_compute(viscous_force_phi) .or. &
sometimes_compute(visc_work) .or. &
sometimes_compute(curl_viscous_force_r) .or. &
sometimes_compute(curl_viscous_force_r_squared) .or. &
sometimes_compute(curl_viscous_force_theta) .or. &
sometimes_compute(curl_viscous_force_theta_squared)) Then
nvf = nvf + 1
vf_p = nvf
Endif

If (compute_quantity(viscous_pforce_r) .or. &
compute_quantity(visc_work_pp)) Then
If (sometimes_compute(viscous_pforce_r) .or. &
sometimes_compute(visc_work_pp) .or. &
sometimes_compute(curl_viscous_pforce_theta) .or. &
sometimes_compute(curl_viscous_pforce_phi)) Then
nvf = nvf+1
vfp_r = nvf
Endif

If (compute_quantity(viscous_pforce_theta) .or. &
compute_quantity(visc_work_pp)) Then
If (sometimes_compute(viscous_pforce_theta) .or. &
sometimes_compute(visc_work_pp) .or. &
sometimes_compute(curl_viscous_pforce_r) .or. &
sometimes_compute(curl_viscous_pforce_phi)) Then
nvf = nvf + 1
vfp_t = nvf
Endif

If (compute_quantity(viscous_pforce_phi) .or. &
compute_quantity(visc_work_pp)) Then
If (sometimes_compute(viscous_pforce_phi) .or. &
sometimes_compute(visc_work_pp) .or. &
sometimes_compute(curl_viscous_pforce_r) .or. &
sometimes_compute(curl_viscous_pforce_theta)) Then
nvf = nvf + 1
vfp_p = nvf
Endif

If (compute_quantity(viscous_mforce_r) .or. &
compute_quantity(visc_work_mm)) Then
If (sometimes_compute(viscous_mforce_r) .or. &
sometimes_compute(visc_work_mm) .or. &
sometimes_compute(curl_viscous_mforce_theta) .or. &
sometimes_compute(curl_viscous_mforce_phi)) Then
nvf = nvf+1
vfm_r = nvf
Endif

If (compute_quantity(viscous_mforce_theta) .or. &
compute_quantity(visc_work_mm)) Then
If (sometimes_compute(viscous_mforce_theta) .or. &
sometimes_compute(visc_work_mm) .or. &
sometimes_compute(curl_viscous_mforce_r) .or. &
sometimes_compute(curl_viscous_mforce_phi)) Then
nvf = nvf + 1
vfm_t = nvf
Endif

If (compute_quantity(viscous_mforce_phi) .or. &
compute_quantity(samom_diffusion) .or. &
compute_quantity(visc_work_mm)) Then
If (sometimes_compute(viscous_mforce_phi) .or. &
sometimes_compute(samom_diffusion) .or. &
sometimes_compute(visc_work_mm) .or. &
sometimes_compute(curl_viscous_mforce_r) .or. &
sometimes_compute(curl_viscous_mforce_theta)) Then
nvf = nvf + 1
vfm_p = nvf
Endif
Expand Down Expand Up @@ -598,6 +622,10 @@ Subroutine Viscous_Force(buffer)

If (compute_quantity(viscous_force_r) .or. &
compute_quantity(visc_work) .or. &
compute_quantity(curl_viscous_force_theta) .or. &
compute_quantity(curl_viscous_force_theta_squared) .or. &
compute_quantity(curl_viscous_force_phi) .or. &
compute_quantity(curl_viscous_force_phi_squared) .or. &
compute_quantity(viscous_mforce_r) ) Then

DO_PSI
Expand Down Expand Up @@ -629,7 +657,11 @@ Subroutine Viscous_Force(buffer)

!Theta-direction; Full
If (compute_quantity(viscous_force_theta) .or. &
compute_quantity(visc_work)) Then
compute_quantity(visc_work) .or. &
compute_quantity(curl_viscous_force_r) .or. &
compute_quantity(curl_viscous_force_r_squared) .or. &
compute_quantity(curl_viscous_force_phi) .or. &
compute_quantity(curl_viscous_force_phi_squared)) Then

DO_PSI
! first, compute all the terms multiplied by mu
Expand All @@ -655,7 +687,11 @@ Subroutine Viscous_Force(buffer)

!Phi-direction
If (compute_quantity(viscous_force_phi) .or. &
compute_quantity(visc_work)) Then
compute_quantity(visc_work) .or. &
compute_quantity(curl_viscous_force_r) .or. &
compute_quantity(curl_viscous_force_r_squared) .or. &
compute_quantity(curl_viscous_force_theta) .or. &
compute_quantity(curl_viscous_force_theta_squared)) Then

DO_PSI
del2u = DDBUFF(PSI,dvpdrdr)+Two_Over_R(r)*buffer(PSI,dvpdr)
Expand Down Expand Up @@ -684,7 +720,9 @@ Subroutine Viscous_Force(buffer)

! r-direction; fluctuating
If (compute_quantity(viscous_pforce_r) .or. &
compute_quantity(visc_work_pp)) Then
compute_quantity(visc_work_pp) .or. &
compute_quantity(curl_viscous_pforce_theta) .or. &
compute_quantity(curl_viscous_pforce_phi)) Then

DO_PSI
! first, compute all the terms multiplied by mu
Expand Down Expand Up @@ -715,7 +753,9 @@ Subroutine Viscous_Force(buffer)

!Theta-direction; Fluctuating
If (compute_quantity(viscous_pforce_theta) .or. &
compute_quantity(visc_work_pp)) Then
compute_quantity(visc_work_pp) .or. &
compute_quantity(curl_viscous_pforce_r) .or. &
compute_quantity(curl_viscous_pforce_phi)) Then

DO_PSI
! first, compute all the terms multiplied by mu
Expand All @@ -741,7 +781,9 @@ Subroutine Viscous_Force(buffer)

!Phi-direction (fluctuating)
If (compute_quantity(viscous_pforce_phi) .or. &
compute_quantity(visc_work_pp)) Then
compute_quantity(visc_work_pp) .or. &
compute_quantity(curl_viscous_pforce_r) .or. &
compute_quantity(curl_viscous_pforce_theta)) Then

DO_PSI
del2u = d2_fbuffer(PSI,dvpdrdr)+Two_Over_R(r)*fbuffer(PSI,dvpdr)
Expand Down Expand Up @@ -770,7 +812,9 @@ Subroutine Viscous_Force(buffer)

! r-direction; mean
If (compute_quantity(viscous_mforce_r) .or. &
compute_quantity(visc_work_mm)) Then
compute_quantity(visc_work_mm) .or. &
compute_quantity(curl_viscous_mforce_theta) .or. &
compute_quantity(curl_viscous_mforce_phi)) Then

DO_PSI
! first, compute all the terms multiplied by mu
Expand Down Expand Up @@ -799,7 +843,9 @@ Subroutine Viscous_Force(buffer)

!Theta-direction; Mean
If (compute_quantity(viscous_mforce_theta) .or. &
compute_quantity(visc_work_mm)) Then
compute_quantity(visc_work_mm) .or. &
compute_quantity(curl_viscous_mforce_r) .or. &
compute_quantity(curl_viscous_mforce_phi)) Then

DO_PSI
! first, compute all the terms multiplied by mu
Expand All @@ -826,7 +872,9 @@ Subroutine Viscous_Force(buffer)
!Phi-direction (mean)
If (compute_quantity(viscous_mforce_phi) .or. &
compute_quantity(samom_diffusion) .or. &
compute_quantity(visc_work_mm)) Then
compute_quantity(visc_work_mm) .or. &
compute_quantity(curl_viscous_mforce_r) .or. &
compute_quantity(curl_viscous_mforce_theta)) Then

DO_PSI
del2u = d2_m0(PSI2,dvpdrdr)+Two_Over_R(r)*m0_values(PSI2,dvpdr)
Expand Down Expand Up @@ -856,11 +904,4 @@ Subroutine Viscous_Force(buffer)

End Subroutine Viscous_Force

Subroutine Finalize_Viscous_Force()
Implicit None
If (nvf .gt. 0) Then
DeAllocate(vforce_buffer)
Endif
End Subroutine Finalize_Viscous_Force

End Module Diagnostics_Linear_Forces
Loading