From 09bc5e312083c9c88c3a9219d18a8227e5921f62 Mon Sep 17 00:00:00 2001 From: Ayesha Sarwar <2831664s@euclid-32.maths.gla.ac.uk> Date: Fri, 17 Jul 2026 03:09:59 +0100 Subject: [PATCH 01/13] Add curl momentum diagnostics --- CONTRIBUTING.md | 1 + src/Diagnostics/Diagnostics_Curl_Momentum.F90 | 471 ++++++++++++++++++ 2 files changed, 472 insertions(+) create mode 100644 src/Diagnostics/Diagnostics_Curl_Momentum.F90 diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 25b3d21c9..d10426565 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -1,3 +1,4 @@ +I am testing the pull request # Contributing to Rayleigh Rayleigh is a community project that lives by the participation of its members — i.e., including you! It is our goal to build an inclusive diff --git a/src/Diagnostics/Diagnostics_Curl_Momentum.F90 b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 new file mode 100644 index 000000000..317e6eb9e --- /dev/null +++ b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 @@ -0,0 +1,471 @@ +Module Diagnostics_Custom + Use Diagnostics_Base + Implicit None +Contains + + Subroutine Custom_MHD_Diagnostics(buffer) + Implicit None + Real*8, Intent(InOut) :: buffer(1:,my_r%min:,my_theta%min:,1:) + Real*8 :: pfactor(my_r%min:my_r%max) + Integer :: r,k, t + pfactor(my_r%min:my_r%max) = ref%dpdr_w_term(my_r%min:my_r%max) & + /ref%density(my_r%min:my_r%max) + + + + !!!!!!!!!!!!!!!!!!!!!!!!! Advection Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + + If (compute_quantity(curl_v_grad_v_r) .or. compute_quantity(curl_v_grad_v_r_squared)) Then + If (compute_quantity(curl_v_grad_v_r)) Then + DO_PSI + qty(PSI) = one_over_r(r) * ref%density(r) * (buffer(PSI,vr) * one_over_r(r) * DDBUFF(PSI,dvpdrdt) + & + buffer(PSI,dvrdt) * one_over_r(r) * buffer(PSI,dvpdr) + buffer(PSI,vtheta) * & + one_over_r(r) * DDBUFF(PSI,dvpdtdt) + buffer(PSI,dvtdt) * one_over_r(r) * & + buffer(PSI,dvpdt) + buffer(PSI,vphi) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvpdtdp) - & + csctheta(t) * csctheta(t) * costheta(t) * buffer(PSI,dvpdp) + buffer(PSI,dvrdt) + & + cottheta(t) * buffer(PSI,dvtdt) - csctheta(t) * & + csctheta(t) * buffer(PSI,vtheta)) + buffer(PSI,dvpdt) * one_over_r(r) * (csctheta(t) * & + buffer(PSI,dvpdp) + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) + & + one_over_r(r) * costheta(t) * csctheta(t) * ref%density(r) * (buffer(PSI,vr) * & + one_over_r(r) * buffer(PSI,dvpdr) + buffer(PSI,vtheta) * one_over_r(r) * & + buffer(PSI,dvpdt) + buffer(PSI,vphi) * one_over_r(r) * (csctheta(t) * & + buffer(PSI,dvpdp) + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) - & + one_over_r(r) * csctheta(t) * ref%density(r) * (buffer(PSI,vr) * DDBUFF(PSI,dvtdrdp) + & + buffer(PSI,dvrdp) * buffer(PSI,dvtdr) + buffer(PSI,vtheta) * one_over_r(r) * & + (DDBUFF(PSI,dvtdtdp) + buffer(PSI,dvrdp)) + buffer(PSI,dvtdp) * one_over_r(r) * & + (buffer(PSI,dvtdt) + buffer(PSI,vr)) + buffer(PSI,vphi) * one_over_r(r) * & + (csctheta(t) * DDBUFF(PSI,dvtdpdp) - cottheta(t) * buffer(PSI,dvpdp)) + & + buffer(PSI,dvpdp) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvtdp) - & + cottheta(t) * buffer(PSI,vphi))) + END_DO + Call Add_Quantity(qty) + Endif + + ! New squared + sqrt quantity (i.e. |curl_v_grad_v_r|) + If (compute_quantity(curl_v_grad_v_r_squared)) Then + DO_PSI + qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + END_DO + Call Add_Quantity(qty) + Endif + + Endif + + + If (compute_quantity(curl_v_grad_v_theta) .or. compute_quantity(curl_v_grad_v_theta_squared)) Then + If (compute_quantity(curl_v_grad_v_theta)) Then + DO_PSI + qty(PSI) = one_over_r(r) * csctheta(t) * ref%density(r) * ( buffer(PSI,vr) * DDBUFF(PSI,dvrdrdp) + & + buffer(PSI,dvrdp) * buffer(PSI,dvrdr) + buffer(PSI,vtheta) * one_over_r(r) * & + (DDBUFF(PSI,dvrdtdp) - buffer(PSI,dvtdp)) + buffer(PSI,dvtdp) * one_over_r(r) * & + (buffer(PSI,dvrdt) - buffer(PSI,vtheta)) + buffer(PSI,vphi) * one_over_r(r) * & + (csctheta(t) * DDBUFF(PSI,dvrdpdp) - buffer(PSI,dvpdp)) + buffer(PSI,dvpdp) * & + one_over_r(r) * ( csctheta(t) * buffer(PSI,dvrdp) - buffer(PSI,vphi))) - & + one_over_r(r) * ref%density(r) * (buffer(PSI,vr) * one_over_r(r) * buffer(PSI,dvpdr) + & + buffer(PSI,vtheta) * one_over_r(r) * buffer(PSI,dvpdt) + buffer(PSI,vphi) * one_over_r(r) * & + (csctheta(t) * buffer(PSI,dvpdp) + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) - & + ref%density(r) * (buffer(PSI,vr) * one_over_r(r) * DDBUFF(PSI,dvpdrdr) + buffer(PSI,dvrdr) * & + one_over_r(r) * buffer(PSI,dvpdr) - buffer(PSI,vr) * one_over_r(r) * one_over_r(r) * & + buffer(PSI,dvpdr) + buffer(PSI,vtheta) * one_over_r(r) * DDBUFF(PSI,dvpdrdt) + & + buffer(PSI,dvtdr) * one_over_r(r) * buffer(PSI,dvpdt) - buffer(PSI,vtheta) * one_over_r(r) * & + one_over_r(r) * buffer(PSI,dvpdt) + & + buffer(PSI,vtheta) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvpdrdp) + & + buffer(PSI,dvrdr) + cottheta(t) * buffer(PSI,dvtdr)) + & + buffer(PSI,dvpdr) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvpdp) + buffer(PSI,vr) + & + cottheta(t) * buffer(PSI,vtheta)) - buffer(PSI,vphi) * one_over_r(r) * one_over_r(r) * & + (csctheta(t) * buffer(PSI,dvpdp) + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) - & + ref%dlnrho(r) * & + ref%density(r) * (buffer(PSI,vr) * one_over_r(r) * buffer(PSI,dvpdr) + & + buffer(PSI,vtheta) * one_over_r(r) * buffer(PSI,dvpdt) + buffer(PSI,vphi) * one_over_r(r) * & + (csctheta(t) * buffer(PSI,dvpdp) + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) + END_DO + Call Add_Quantity(qty) + Endif + If (compute_quantity(curl_v_grad_v_theta_squared)) Then + DO_PSI + qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + END_DO + Call Add_Quantity(qty) + Endif + + Endif + + If (compute_quantity(curl_v_grad_v_phi) .or. compute_quantity(curl_v_grad_v_phi_squared)) Then + If (compute_quantity(curl_v_grad_v_theta)) Then + DO_PSI + qty(PSI) = ref%density(r) * (buffer(PSI,vr) * DDBUFF(PSI,dvtdrdr) + buffer(PSI,dvrdr) * buffer(PSI,dvtdr) + & + buffer(PSI,vtheta) * one_over_r(r) * (DDBUFF(PSI,dvtdrdt) + buffer(PSI,dvrdr)) + & + buffer(PSI,dvtdr) * one_over_r(r) * (buffer(PSI,dvtdt) + buffer(PSI,vr)) - & + buffer(PSI,vtheta) * one_over_r(r) * one_over_r(r) * (buffer(PSI,dvtdt) + buffer(PSI,vr)) + & + buffer(PSI,vtheta) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvtdrdp) - cottheta(t) * & + buffer(PSI,dvpdr)) + buffer(PSI,dvpdr) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvtdp) - & + cottheta(t) * buffer(PSI,vphi)) - buffer(PSI,vphi) * one_over_r(r) * one_over_r(r) * & + (csctheta(t) * buffer(PSI,dvtdp) - cottheta(t) * buffer(PSI,vphi))) + & + ref%dlnrho(r) * ref%density(r) * (buffer(PSI,vr) * buffer(PSI,dvtdr) + buffer(PSI,vtheta) * & + one_over_r(r) * (buffer(PSI,dvtdt) + buffer(PSI,vr)) + buffer(PSI,vphi) * one_over_r(r) * & + (csctheta(t) * buffer(PSI,dvtdp) - cottheta(t) * buffer(PSI,vphi))) + & + one_over_r(r) * ref%density(r) * (buffer(PSI,vr) * buffer(PSI,dvtdr) + buffer(PSI,vtheta) * & + one_over_r(r) * (buffer(PSI,dvtdt) + buffer(PSI,vr)) + buffer(PSI,vphi) * one_over_r(r) * & + (csctheta(t) * buffer(PSI,dvtdp) - cottheta(t) * buffer(PSI,vphi))) - & + one_over_r(r) * ref%density(r) * (buffer(PSI,vr) * DDBUFF(PSI,dvrdrdt) + buffer(PSI,dvrdt) * & + buffer(PSI,dvrdr) + buffer(PSI,vtheta) * one_over_r(r) * (DDBUFF(PSI,dvrdtdt) + buffer(PSI,dvtdt)) + & + buffer(PSI,dvtdt) * one_over_r(r) * (buffer(PSI,dvrdt) - buffer(PSI,vtheta)) + & + buffer(PSI,vphi) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvrdtdp) - csctheta(t) * csctheta(t) * & + costheta(t) * buffer(PSI,dvrdp) - buffer(PSI,dvpdt)) + buffer(PSI,dvpdt) * one_over_r(r) * & + (csctheta(t) * buffer(PSI,dvrdp) - buffer(PSI,vphi))) + + END_DO + Call Add_Quantity(qty) + Endif + If (compute_quantity(curl_v_grad_v_phi_squared)) Then + DO_PSI + qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + END_DO + Call Add_Quantity(qty) + Endif + + Endif + + !!!!!!!!!!!!!!!!!!!!!!!!!! Buoyancy Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + If (compute_quantity(curl_buoyancy_force_theta) .or. compute_quantity(curl_buoyancy_force_theta_squared)) Then + If (compute_quantity(curl_buoyancy_force_theta)) Then + DO_PSI + qty(PSI) = ref%Buoyancy_Coeff(r) * (csctheta(t) * & + radius(r) * buffer(PSI,dtdp)) ! since dtdp = (1/r)*d(temperature or entropy)/dphi + END_DO + Call Add_Quantity(qty) + Endif + + ! New squared + sqrt quantity (i.e. |curl_buoyancy_force_theta|) + If (compute_quantity(curl_buoyancy_force_theta_squared)) Then + DO_PSI + qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + END_DO + Call Add_Quantity(qty) + Endif + + Endif + + If (compute_quantity(curl_buoyancy_force_phi) .or. compute_quantity(curl_buoyancy_force_phi_squared)) Then + If (compute_quantity(curl_buoyancy_force_phi)) Then + DO_PSI + qty(PSI) = -ref%Buoyancy_Coeff(r) * ( radius(r) * buffer(PSI,dtdt)) ! since dtdt = (1/r)*d(temperature or $ + END_DO + Call Add_Quantity(qty) + Endif + + ! New squared + sqrt quantity (i.e. |curl_buoyancy_force_phi|) + If (compute_quantity(curl_buoyancy_force_phi_squared)) Then + DO_PSI + qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + END_DO + Call Add_Quantity(qty) + Endif + + Endif + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Magnetic Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + If (compute_quantity(curl_j_cross_b_r) .or. compute_quantity(curl_j_cross_b_r_squared)) Then + If (compute_quantity(curl_j_cross_b_r)) Then + DO_PSI + qty(PSI) = ref%Lorentz_Coeff*(one_over_r(r) * one_over_r(r) * buffer(PSI,dbtdt) * buffer(PSI,dbpdt) + & + one_over_r(r) * one_over_r(r) * buffer(PSI,btheta) * DDBUFF(PSI,dbpdtdt) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * buffer(PSI,btheta) * & + buffer(PSI,bphi) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,dbtdt) * & + buffer(PSI,bphi) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,btheta) * & + buffer(PSI,dbpdt) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * buffer(PSI,dbtdt) * buffer(PSI,dbtdp) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * buffer(PSI,btheta) * DDBUFF(PSI,dbtdtdp) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * buffer(PSI,btheta) * & + buffer(PSI,dbtdp) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * & + buffer(PSI,br) * buffer(PSI,dbrdp) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,dbrdt) * buffer(PSI,dbrdp) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,br) * DDBUFF(PSI,dbrdtdp) + & + one_over_r(r) * one_over_r(r) * buffer(PSI,dbrdt) * buffer(PSI,bphi) + & + one_over_r(r) * one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbpdt) + & + one_over_r(r) * buffer(PSI,dbrdt) * buffer(PSI,dbpdr) + & + one_over_r(r) * buffer(PSI,br) * DDBUFF(PSI,dbpdrdt) + & !(first part ended) + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & + buffer(PSI,btheta) * buffer(PSI,dbpdt) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * costheta(t) * & + buffer(PSI,btheta) * buffer(PSI,bphi) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * & + buffer(PSI,btheta) * buffer(PSI,dbtdp) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * & + buffer(PSI,br) * buffer(PSI,dbrdp) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & + buffer(PSI,br) * buffer(PSI,bphi) + & + one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,br) * buffer(PSI,dbpdr) - & !(second part ended) + one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbtdp) - & + one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,dbrdp) - & + buffer(PSI,dbrdp) * buffer(PSI,dbtdr) - & + buffer(PSI,br) * DDBUFF(PSI,dbtdrdp) + & + buffer(PSI,dbrdp) * one_over_r(r) * buffer(PSI,dbrdt) + & + buffer(PSI,br) * one_over_r(r) * DDBUFF(PSI,dbrdtdp) + & + one_over_r(r) * buffer(PSI,dbpdp) * buffer(PSI,dbpdt) + & + one_over_r(r) * buffer(PSI,bphi) * DDBUFF(PSI,dbpdtdp) + & + 2 * one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,bphi) * buffer(PSI,dbpdp) - & + buffer(PSI,dbpdp) * buffer(PSI,dbtdp) - & + buffer(PSI,bphi) * DDBUFF(PSI,dbtdpdp)) + END_DO + Call Add_Quantity(qty) + Endif + + ! New squared + sqrt quantity (i.e. |curl_j_cross_b_r|) + If (compute_quantity(curl_j_cross_b_r_squared)) Then + DO_PSI + qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + END_DO + Call Add_Quantity(qty) + Endif + + Endif + + + If (compute_quantity(curl_j_cross_b_theta) .or. compute_quantity(curl_j_cross_b_theta_squared)) Then + If (compute_quantity(curl_j_cross_b_theta)) Then + DO_PSI + qty(PSI) = ref%Lorentz_Coeff*(one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * & + buffer(PSI,dbpdp) * buffer(PSI,dbrdp) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * & + buffer(PSI,bphi) * DDBUFF(PSI,dbrdpdp) - & + 2 * one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,bphi) * buffer(PSI,dbpdp) - & + one_over_r(r) * csctheta(t) * buffer(PSI,dbpdp) * buffer(PSI,dbpdr) - & + one_over_r(r) * csctheta(t) * buffer(PSI,bphi) * DDBUFF(PSI,dbpdrdp) - & + 2 * one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,btheta) * buffer(PSI,dbtdp) - & + one_over_r(r) * csctheta(t) * buffer(PSI,dbtdp) * buffer(PSI,dbtdr) - & + one_over_r(r) * csctheta(t) * buffer(PSI,btheta) * DDBUFF(PSI,dbtdrdp) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,dbtdp) * buffer(PSI,dbrdt) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,btheta) * DDBUFF(PSI,dbrdtdp) - & !(first part ended) + one_over_r(r) * one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,dbpdt) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & + buffer(PSI,btheta) * buffer(PSI,bphi) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,btheta) * buffer(PSI,dbtdp) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * & + buffer(PSI,br) * buffer(PSI,dbrdp) - & + one_over_r(r) * one_over_r(r) * buffer(PSI,br) * buffer(PSI,bphi) - & + one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbpdr) + & !(second part ended) + one_over_r(r) * one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,dbpdt) - & + one_over_r(r) * buffer(PSI,dbtdr) * buffer(PSI,dbpdt) - & + one_over_r(r) * buffer(PSI,btheta) * DDBUFF(PSI,dbpdrdt) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & + buffer(PSI,btheta) * buffer(PSI,bphi) - & + one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,dbtdr) * buffer(PSI,bphi) - & + one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,btheta) * buffer(PSI,dbpdr) + & + one_over_r(r) * csctheta(t) * buffer(PSI,dbtdr) * buffer(PSI,dbtdp) - & + buffer(PSI,btheta) * one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,dbtdr) + & + one_over_r(r) * csctheta(t) * buffer(PSI,btheta) * DDBUFF(PSI,dbtdrdp) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,br) * buffer(PSI,dbrdp) + & + one_over_r(r) * csctheta(t) * buffer(PSI,dbrdr) * buffer(PSI,dbrdp) + & + one_over_r(r) * csctheta(t) * buffer(PSI,br) * DDBUFF(PSI,dbrdrdp) + & + one_over_r(r) * one_over_r(r) * buffer(PSI,br) * buffer(PSI,bphi) - & + one_over_r(r) * buffer(PSI,dbrdr) * buffer(PSI,bphi) - & + one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbpdr) - & + buffer(PSI,dbrdr) * buffer(PSI,dbpdr) - & + buffer(PSI,br) * DDBUFF(PSI,dbpdrdr)) + END_DO + Call Add_Quantity(qty) + Endif + + ! New squared + sqrt quantity (i.e. |curl_j_cross_b_theta|) + If (compute_quantity(curl_j_cross_b_theta_squared)) Then + DO_PSI + qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + END_DO + Call Add_Quantity(qty) + Endif + + Endif + + + If (compute_quantity(curl_j_cross_b_phi) .or. compute_quantity(curl_j_cross_b_phi_squared)) Then + If (compute_quantity(curl_j_cross_b_phi)) Then + DO_PSI + qty(PSI) = ref%Lorentz_Coeff*(one_over_r(r) * one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,br) + & + one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbtdr) - & + one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbrdt) - & + one_over_r(r) * one_over_r(r) * buffer(PSI,bphi) * buffer(PSI,dbpdt) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * (buffer(PSI, bphi))**2 + & + one_over_r(r) * buffer(PSI,bphi) * buffer(PSI,dbtdp) - & !(first part ended) + one_over_r(r) * one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,br) + & + one_over_r(r) * buffer(PSI,dbtdr) * buffer(PSI,br) + & + one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,dbrdr) + & + buffer(PSI,dbrdr) * buffer(PSI,dbtdr) + & + buffer(PSI,br) * DDBUFF(PSI,dbtdrdr) - & + buffer(PSI,dbrdr) * buffer(PSI,dbrdt) - & + one_over_r(r) * buffer(PSI,br) * DDBUFF(PSI,dbrdrdt) + & + one_over_r(r) * one_over_r(r) * buffer(PSI,dbrdt) + & + one_over_r(r) * one_over_r(r) * buffer(PSI,bphi) * buffer(PSI,dbpdt) - & + one_over_r(r) * buffer(PSI,dbpdr) * buffer(PSI,dbpdt) - & + one_over_r(r) * buffer(PSI,bphi) * DDBUFF(PSI,dbpdrdt) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & + (buffer(PSI, bphi))**2 - & + 2 * one_over_r(r) * csctheta(t) * costheta(t) * & + buffer(PSI, bphi) * buffer(PSI,dbpdr) + & + one_over_r(r) * csctheta(t) * buffer(PSI, dbpdr) * buffer(PSI, dbtdp) - & + buffer(PSI, bphi) * one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI, dbtdp) + & + one_over_r(r) * csctheta(t) * buffer(PSI, bphi) * DDBUFF(PSI,dbtdrdp) + & !(second part ended) + one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * & + buffer(PSI, bphi) * buffer(PSI,dbrdp) - & + one_over_r(r) * csctheta(t) * buffer(PSI,dbpdt) * buffer(PSI,dbrdp) - & + one_over_r(r) * csctheta(t) * buffer(PSI,bphi) * DDBUFF(PSI,dbrdtdp) + & + 2 * one_over_r(r) * buffer(PSI,bphi) * buffer(PSI,dbpdt) + & + buffer(PSI,dbpdt) * buffer(PSI,dbpdr) + & + buffer(PSI,bphi) * DDBUFF(PSI,dbpdrdt) + & + 2 * one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,dbtdt) + & + buffer(PSI,dbtdt) * buffer(PSI,dbtdr) + & + buffer(PSI,btheta) * DDBUFF(PSI,dbtdrdt) - & + one_over_r(r) * buffer(PSI,dbtdt) * buffer(PSI,dbrdt) - & + one_over_r(r) * buffer(PSI,btheta) * DDBUFF(PSI,dbrdtdt)) + END_DO + Call Add_Quantity(qty) + Endif + + ! New squared + sqrt quantity (i.e. |curl_j_cross_b_phi|) + If (compute_quantity(curl_j_cross_b_phi_squared)) Then + DO_PSI + qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + END_DO + Call Add_Quantity(qty) + Endif + + Endif + + !!!!!!!!!!!!!!!!!!!!!!!!!!!! Coriolis Force !!!!!!!!!!!!!!!!!!!!!!!!!! + + If (compute_quantity(curl_coriolis_force_r) .or. compute_quantity(curl_coriolis_force_r_squared)) Then + If (compute_quantity(curl_coriolis_force_r)) Then + DO_PSI + qty(PSI) = - ref%Coriolis_Coeff * ref%density(r) * one_over_r(r) * & + (-sintheta(t) * buffer(PSI,vtheta) + cottheta(t) * csctheta(t) * buffer(PSI,vtheta) + & + costheta(t) * buffer(PSI,dvtdt) + 2 * costheta(t) * buffer(PSI,vr) + sintheta(t) * & + buffer(PSI,dvrdt) + cottheta(t) * buffer(PSI,dvpdt)) + END_DO + Call Add_Quantity(qty) + Endif + + ! New squared + sqrt quantity (i.e. |curl_coriolis_force_r|) + If (compute_quantity(curl_coriolis_force_r_squared)) Then + DO_PSI + qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + END_DO + Call Add_Quantity(qty) + Endif + + Endif + + If (compute_quantity(curl_coriolis_force_theta) .or. compute_quantity(curl_coriolis_force_theta_squared)) Then + If (compute_quantity(curl_coriolis_force_theta)) Then + DO_PSI + qty(PSI) = ref%Coriolis_Coeff * ref%density(r) * (one_over_r(r) * (buffer(PSI,dvpdp) + & + costheta(t) * buffer(PSI,vtheta) + sintheta(t) * buffer(PSI,vr)) + costheta(t) * & + buffer(PSI,dvtdr) + sintheta(t) * buffer(PSI,dvrdr) + ref%dlnrho(r) * costheta(t) * & + buffer(PSI,vtheta) + ref%dlnrho(r) * sintheta(t) * buffer(PSI,vr)) + END_DO + Call Add_Quantity(qty) + Endif + + ! New squared + sqrt quantity (i.e. |curl_coriolis_force_theta|) + If (compute_quantity(curl_coriolis_force_theta_squared)) Then + DO_PSI + qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + END_DO + Call Add_Quantity(qty) + Endif + + Endif + + If (compute_quantity(curl_coriolis_force_phi) .or. compute_quantity(curl_coriolis_force_phi)) Then + If (compute_quantity(curl_coriolis_force_phi)) Then + DO_PSI + qty(PSI) = ref%Coriolis_Coeff * ref%density(r) * (ref%dlnrho(r) * costheta(t) * buffer(PSI,vphi) + & + costheta(t) * buffer(PSI,dvpdr) - one_over_r(r) * sintheta(t) * buffer(PSI,dvpdr)) + + END_DO + Call Add_Quantity(qty) + Endif + + ! New squared + sqrt quantity (i.e. |curl_coriolis_force_phi|) + If (compute_quantity(curl_coriolis_force_phi_squared)) Then + DO_PSI + qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + END_DO + Call Add_Quantity(qty) + Endif + + Endif + + + + !!!!!!!!!!!!!!!!!!!!!!!!!!! Pressure Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + + If (compute_quantity(curl_pressure_force_r) .or. compute_quantity(curl_pressure_force_r_squared)) Then + If (compute_quantity(curl_pressure_force_r)) Then + DO_PSI + qty(PSI) = pfactor(r) * & + OneOverRSquared(r) * csctheta(t) * (DDBUFF(PSI,dvpdtdp)-DDBUFF(PSI,dvpdtdp)) + END_DO + Call Add_Quantity(qty) + Endif + + ! New squared + sqrt quantity (i.e. |curl_viscous_force_r|) + If (compute_quantity(curl_pressure_force_r_squared)) Then + DO_PSI + qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + END_DO + Call Add_Quantity(qty) + Endif + Endif + + If (compute_quantity(curl_pressure_force_theta) .or. compute_quantity(curl_pressure_force_theta_squared)) Then + If (compute_quantity(curl_pressure_force_theta)) Then + DO_PSI + qty(PSI) = pfactor(r) * & + OneOverRSquared(r) * csctheta(t) * (DDBUFF(PSI,dvpdrdp) + & + ref%dlnrho(r) * buffer(PSI,dpdp) + DDBUFF(PSI,dvpdrdp)) + END_DO + Call Add_Quantity(qty) + Endif + + ! New squared + sqrt quantity (i.e. |curl_pressure_force_theta|) + If (compute_quantity(curl_pressure_force_theta_squared)) Then + DO_PSI + qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + END_DO + Call Add_Quantity(qty) + Endif + Endif + + + + If (compute_quantity(curl_pressure_force_phi) .or. compute_quantity(curl_pressure_force_phi_squared)) Then + If (compute_quantity(curl_pressure_force_phi)) Then + DO_PSI + qty(PSI) = pfactor(r) * & + one_over_r(r) * csctheta(t) * (-DDBUFF(PSI,dvpdrdt) + DDBUFF(PSI,dvpdrdt) - & + ref%dlnrho(r) * buffer(PSI,dpdt)) + END_DO + Call Add_Quantity(qty) + Endif + + ! New squared + sqrt quantity (i.e. |curl_viscous_force_phi|) + If (compute_quantity(curl_pressure_force_phi_squared)) Then + DO_PSI + qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + END_DO + Call Add_Quantity(qty) + Endif + Endif + + \ No newline at end of file From 19c9e4e95992b7d87bc5c16a60356d1bf9161470 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Thu, 16 Jul 2026 20:23:58 -0600 Subject: [PATCH 02/13] Removing line --- CONTRIBUTING.md | 1 - 1 file changed, 1 deletion(-) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index d10426565..25b3d21c9 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -1,4 +1,3 @@ -I am testing the pull request # Contributing to Rayleigh Rayleigh is a community project that lives by the participation of its members — i.e., including you! It is our goal to build an inclusive From 7db7fede043b35ef7eb33a03b2a060d39733d4cc Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Fri, 17 Jul 2026 09:15:50 -0600 Subject: [PATCH 03/13] Adding a placeholder file to contain the curl of the momentum equation quantity codes (currently is mostly a copy of the momentum equation codes so not valid). --- .../curl_momentum_equation_codes.F | 123 ++++++++++++++++++ 1 file changed, 123 insertions(+) create mode 100644 src/Diagnostics/curl_momentum_equation_codes.F diff --git a/src/Diagnostics/curl_momentum_equation_codes.F b/src/Diagnostics/curl_momentum_equation_codes.F new file mode 100644 index 000000000..b904b413b --- /dev/null +++ b/src/Diagnostics/curl_momentum_equation_codes.F @@ -0,0 +1,123 @@ +! +! Copyright (C) 2018 by the authors of the RAYLEIGH code. +! +! This file is part of RAYLEIGH. +! +! RAYLEIGH is free software; you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation; either version 3, or (at your option) +! any later version. +! +! RAYLEIGH is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with RAYLEIGH; see the file LICENSE. If not see +! . +! + + !//////////////////////// Advection Terms //////////////////// + ! Reynolds decomposition about the azimuthal mean may also be output + ! NOTE: ADVECTION TERMS ARE SCALED BY DENSITY (so that they represent a force density) + + Integer, Parameter :: curl_mom_eq_off = mom_eq_off+100 ! :OFFSET CODE: + Integer, Parameter :: v_grad_v_r = curl_mom_eq_off+1 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_r$ + Integer, Parameter :: v_grad_v_theta = curl_mom_eq_off+2 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\theta$ + Integer, Parameter :: v_grad_v_phi = curl_mom_eq_off+3 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\phi$ + + Integer, Parameter :: vp_grad_vm_r = curl_mom_eq_off+4 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_r$ + Integer, Parameter :: vp_grad_vm_theta = curl_mom_eq_off+5 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\theta$ + Integer, Parameter :: vp_grad_vm_phi = curl_mom_eq_off+6 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\phi$ + + Integer, Parameter :: vm_grad_vp_r = curl_mom_eq_off+7 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_r$ + Integer, Parameter :: vm_grad_vp_theta = curl_mom_eq_off+8 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\theta$ + Integer, Parameter :: vm_grad_vp_phi = curl_mom_eq_off+9 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\phi$ + + Integer, Parameter :: vp_grad_vp_r = curl_mom_eq_off+10 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_r$ + Integer, Parameter :: vp_grad_vp_theta = curl_mom_eq_off+11 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\theta$ + Integer, Parameter :: vp_grad_vp_phi = curl_mom_eq_off+12 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\phi$ + + Integer, Parameter :: vm_grad_vm_r = curl_mom_eq_off+13 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_r$ + Integer, Parameter :: vm_grad_vm_theta = curl_mom_eq_off+14 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\theta$ + Integer, Parameter :: vm_grad_vm_phi = curl_mom_eq_off+15 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\phi$ + + !///////////////////////////////////////////////////////////// + ! Linear Forces + ! Note: the pressure gradient diagnostic codes are above + + Integer, Parameter :: buoyancy_force = curl_mom_eq_off+16 ! :tex: $c_2\mathrm{f}_2\Theta$ + Integer, Parameter :: buoyancy_pforce = curl_mom_eq_off+17 ! :tex: $c_2\mathrm{f}_2\Theta'$ + Integer, Parameter :: buoyancy_mforce = curl_mom_eq_off+18 ! :tex: $c_2\mathrm{f}_2\overline{\Theta}$ + + Integer, Parameter :: Coriolis_Force_r = curl_mom_eq_off+19 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\boldsymbol{v}\right]_r$ + Integer, Parameter :: Coriolis_Force_theta = curl_mom_eq_off+20 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\boldsymbol{v}\right]_\theta$ + Integer, Parameter :: Coriolis_Force_phi = curl_mom_eq_off+21 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\boldsymbol{v}\right]_\phi$ + + Integer, Parameter :: Coriolis_pForce_r = curl_mom_eq_off+22 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\boldsymbol{v'}\right]_r$ + Integer, Parameter :: Coriolis_pForce_theta = curl_mom_eq_off+23 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\boldsymbol{v'}\right]_\theta$ + Integer, Parameter :: Coriolis_pForce_phi = curl_mom_eq_off+24 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\boldsymbol{v'}\right]_\phi$ + + Integer, Parameter :: Coriolis_mForce_r = curl_mom_eq_off+25 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\overline{\boldsymbol{v}}\right]_r$ + Integer, Parameter :: Coriolis_mForce_theta = curl_mom_eq_off+26 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\overline{\boldsymbol{v}}\right]_\theta$ + Integer, Parameter :: Coriolis_mForce_phi = curl_mom_eq_off+27 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\overline{\boldsymbol{v}}\right]_\phi$ + + ! Viscous forces + Integer, Parameter :: viscous_Force_r = curl_mom_eq_off+28 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D}}\right]_r$ + Integer, Parameter :: viscous_Force_theta = curl_mom_eq_off+29 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D}}\right]_\theta$ + Integer, Parameter :: viscous_Force_phi = curl_mom_eq_off+30 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D}}\right]_\phi$ + + Integer, Parameter :: viscous_pForce_r = curl_mom_eq_off+31 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D'}}\right]_r$ + Integer, Parameter :: viscous_pForce_theta = curl_mom_eq_off+32 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D'}}\right]_\theta$ + Integer, Parameter :: viscous_pForce_phi = curl_mom_eq_off+33 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D'}}\right]_\phi$ + + Integer, Parameter :: viscous_mForce_r = curl_mom_eq_off+34 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal{D}}}\right]_r$ + Integer, Parameter :: viscous_mForce_theta = curl_mom_eq_off+35 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal{D}}}\right]_\theta$ + Integer, Parameter :: viscous_mForce_phi = curl_mom_eq_off+36 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal{D}}}\right]_\phi$ + + ! Pressure forces + Integer, Parameter :: pressure_Force_r = curl_mom_eq_off+37 ! :tex: $ -c_3\mathrm{f}_1\frac{\partial}{\partial r}\left(\frac{P}{\mathrm{f}_1} \right)$ + Integer, Parameter :: pressure_Force_theta = curl_mom_eq_off+38 ! :tex: $ -c_3\frac{1}{r}\frac{\partial P}{\partial \theta}$ + Integer, Parameter :: pressure_Force_phi = curl_mom_eq_off+39 ! :tex: $ -c_3\frac{1}{r\mathrm{sin}\theta}\frac{\partial P}{\partial \phi}$ + + Integer, Parameter :: pressure_pForce_r = curl_mom_eq_off+40 ! :tex: $ -c_3\mathrm{f}_1\frac{\partial}{\partial r}\left(\frac{P'}{\mathrm{f}_1} \right)$ + Integer, Parameter :: pressure_pForce_theta = curl_mom_eq_off+41 ! :tex: $ -c_3\frac{1}{r}\frac{\partial P'}{\partial \theta}$ + Integer, Parameter :: pressure_pForce_phi = curl_mom_eq_off+42 ! :tex: $ -c_3\frac{1}{r\mathrm{sin}\theta}\frac{\partial P'}{\partial \phi}$ + + Integer, Parameter :: pressure_mForce_r = curl_mom_eq_off+43 ! :tex: $ -c_3\mathrm{f}_1\frac{\partial}{\partial r}\left(\frac{\overline{P}}{\mathrm{f}_1} \right)$ + Integer, Parameter :: pressure_mForce_theta = curl_mom_eq_off+44 ! :tex: $ -c_3\frac{1}{r}\frac{\partial \overline{P}}{\partial \theta}$ + Integer, Parameter :: pressure_mForce_phi = curl_mom_eq_off+45 ! :tex: $ -c_3\frac{1}{r\mathrm{sin}\theta}\frac{\partial \overline{P}}{\partial \phi}$ + + ! ell=0 pressure and buoyancy forces (r-direction only) + ! These are substracted out from the radial terms above + ! because they do not contribute to the dynamics. + ! The ell=0 buoyancy force is determined by the mean entropy profile. + ! The ell=0 pressure force is DEFINED such that it exactly balances the ell=0 buoyancy + Integer, Parameter :: buoyancy_force_ell0 = curl_mom_eq_off+46 ! :tex: $c_2\mathrm{f}_2\Theta_{00}$ + Integer, Parameter :: pressure_force_ell0_r = curl_mom_eq_off+47 ! :tex: $-c_3\mathrm{f}_1\frac{\partial}{\partial r}\left(\frac{P_{00}}{\mathrm{f}_1} \right)$ + + !/////////////////////////// Lorentz Forces /////////////////////////////// + ! ref%Lorentz_Coeff * (del x B) x B + ! ref%Lorentz_Coeff = 1/4pi when dimensional, Pr/(Pr_m E) when nondimesional + ! j (below) is shorthand for ref%Lorentz_Coeff*delxB (not quite the current density) + + Integer, Parameter :: j_cross_b_r = curl_mom_eq_off+48 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B}\right)\times\boldsymbol{B}\right]_r$ + Integer, Parameter :: j_cross_b_theta = curl_mom_eq_off+49 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B}\right)\times\boldsymbol{B}\right]_\theta$ + Integer, Parameter :: j_cross_b_phi = curl_mom_eq_off+50 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B}\right)\times\boldsymbol{B}\right]_\phi$ + + Integer, Parameter :: jp_cross_bm_r = curl_mom_eq_off+51 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B'}\right)\times\overline{\boldsymbol{B}}\right]_r$ + Integer, Parameter :: jp_cross_bm_theta = curl_mom_eq_off+52 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B'}\right)\times\overline{\boldsymbol{B}}\right]_\theta$ + Integer, Parameter :: jp_cross_bm_phi = curl_mom_eq_off+53 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B'}\right)\times\overline{\boldsymbol{B}}\right]_\phi$ + + Integer, Parameter :: jm_cross_bp_r = curl_mom_eq_off+54 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\overline{\boldsymbol{B}}\right)\times\boldsymbol{B'}\right]_r$ + Integer, Parameter :: jm_cross_bp_theta = curl_mom_eq_off+55 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\overline{\boldsymbol{B}}\right)\times\boldsymbol{B'}\right]_\theta$ + Integer, Parameter :: jm_cross_bp_phi = curl_mom_eq_off+56 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\overline{\boldsymbol{B}}\right)\times\boldsymbol{B'}\right]_\phi$ + + Integer, Parameter :: jm_cross_bm_r = curl_mom_eq_off+57 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\overline{\boldsymbol{B}}\right)\times\overline{\boldsymbol{B}}\right]_r$ + Integer, Parameter :: jm_cross_bm_theta = curl_mom_eq_off+58 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\overline{\boldsymbol{B}}\right)\times\overline{\boldsymbol{B}}\right]_\theta$ + Integer, Parameter :: jm_cross_bm_phi = curl_mom_eq_off+59 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\overline{\boldsymbol{B}}\right)\times\overline{\boldsymbol{B}}\right]_\phi$ + + Integer, Parameter :: jp_cross_bp_r = curl_mom_eq_off+60 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B'}\right)\times\boldsymbol{B'}\right]_r$ + Integer, Parameter :: jp_cross_bp_theta = curl_mom_eq_off+61 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B'}\right)\times\boldsymbol{B'}\right]_\theta$ + Integer, Parameter :: jp_cross_bp_phi = curl_mom_eq_off+62 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B'}\right)\times\boldsymbol{B'}\right]_\phi$ From 6c85b5fac7653c80163bf026c217bcb3eb654b7c Mon Sep 17 00:00:00 2001 From: Ayesha Sarwar <2831664s@euclid-32.maths.gla.ac.uk> Date: Fri, 17 Jul 2026 19:12:49 +0100 Subject: [PATCH 04/13] Update curl momentum equation diagnostics --- .../curl_momentum_equation_codes.F | 167 ++++++++++-------- 1 file changed, 93 insertions(+), 74 deletions(-) diff --git a/src/Diagnostics/curl_momentum_equation_codes.F b/src/Diagnostics/curl_momentum_equation_codes.F index b904b413b..b6180face 100644 --- a/src/Diagnostics/curl_momentum_equation_codes.F +++ b/src/Diagnostics/curl_momentum_equation_codes.F @@ -23,101 +23,120 @@ ! NOTE: ADVECTION TERMS ARE SCALED BY DENSITY (so that they represent a force density) Integer, Parameter :: curl_mom_eq_off = mom_eq_off+100 ! :OFFSET CODE: - Integer, Parameter :: v_grad_v_r = curl_mom_eq_off+1 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_r$ - Integer, Parameter :: v_grad_v_theta = curl_mom_eq_off+2 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\theta$ - Integer, Parameter :: v_grad_v_phi = curl_mom_eq_off+3 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\phi$ + Integer, Parameter :: curl_v_grad_v_r = curl_mom_eq_off+1 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_r$ + Integer, Parameter :: curl_v_grad_v_theta = curl_mom_eq_off+2 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\theta$ + Integer, Parameter :: curl_v_grad_v_phi = curl_mom_eq_off+3 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\phi$ - Integer, Parameter :: vp_grad_vm_r = curl_mom_eq_off+4 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_r$ - Integer, Parameter :: vp_grad_vm_theta = curl_mom_eq_off+5 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\theta$ - Integer, Parameter :: vp_grad_vm_phi = curl_mom_eq_off+6 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\phi$ + Integer, Parameter :: curl_v_grad_v_r_squared = curl_mom_eq_off+4 ! :tex: $\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_r\right)^2$ + Integer, Parameter :: curl_v_grad_v_theta_squared = curl_mom_eq_off+5 ! :tex: $\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\theta\right)^2$ + Integer, Parameter :: curl_v_grad_v_phi_squared = curl_mom_eq_off+6 ! :tex: $\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\phi\right)^2$ - Integer, Parameter :: vm_grad_vp_r = curl_mom_eq_off+7 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_r$ - Integer, Parameter :: vm_grad_vp_theta = curl_mom_eq_off+8 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\theta$ - Integer, Parameter :: vm_grad_vp_phi = curl_mom_eq_off+9 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\phi$ +! Integer, Parameter :: curl_vp_grad_vm_r = curl_mom_eq_off+7 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_r$ +! Integer, Parameter :: curl_vp_grad_vm_theta = curl_mom_eq_off+8 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\theta$ +! Integer, Parameter :: curl_vp_grad_vm_phi = curl_mom_eq_off+9 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\phi$ - Integer, Parameter :: vp_grad_vp_r = curl_mom_eq_off+10 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_r$ - Integer, Parameter :: vp_grad_vp_theta = curl_mom_eq_off+11 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\theta$ - Integer, Parameter :: vp_grad_vp_phi = curl_mom_eq_off+12 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\phi$ +! Integer, Parameter :: curl_vm_grad_vp_r = curl_mom_eq_off+10 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_r$ +! Integer, Parameter :: curl_vm_grad_vp_theta = curl_mom_eq_off+11 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\theta$ +! Integer, Parameter :: curl_vm_grad_vp_phi = curl_mom_eq_off+12 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\phi$ - Integer, Parameter :: vm_grad_vm_r = curl_mom_eq_off+13 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_r$ - Integer, Parameter :: vm_grad_vm_theta = curl_mom_eq_off+14 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\theta$ - Integer, Parameter :: vm_grad_vm_phi = curl_mom_eq_off+15 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\phi$ +! Integer, Parameter :: curl_vp_grad_vp_r = curl_mom_eq_off+13 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_r$ +! Integer, Parameter :: curl_vp_grad_vp_theta = curl_mom_eq_off+14 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\theta$ +! Integer, Parameter :: curl_vp_grad_vp_phi = curl_mom_eq_off+15 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\phi$ + +! Integer, Parameter :: curl_vm_grad_vm_r = curl_mom_eq_off+16 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_r$ +! Integer, Parameter :: curl_vm_grad_vm_theta = curl_mom_eq_off+17 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\theta$ +! Integer, Parameter :: curl_vm_grad_vm_phi = curl_mom_eq_off+18 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\phi$ !///////////////////////////////////////////////////////////// - ! Linear Forces + ! Linear forces ! Note: the pressure gradient diagnostic codes are above - Integer, Parameter :: buoyancy_force = curl_mom_eq_off+16 ! :tex: $c_2\mathrm{f}_2\Theta$ - Integer, Parameter :: buoyancy_pforce = curl_mom_eq_off+17 ! :tex: $c_2\mathrm{f}_2\Theta'$ - Integer, Parameter :: buoyancy_mforce = curl_mom_eq_off+18 ! :tex: $c_2\mathrm{f}_2\overline{\Theta}$ + Integer, Parameter :: curl_buoyancy_force_theta = curl_mom_eq_off+19 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\Theta\right)\right]_\theta$ + Integer, Parameter :: curl_buoyancy_force_phi = curl_mom_eq_off+20 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\Theta\right)\right]_\phi$ + + Integer, Parameter :: curl_buoyancy_force_theta_squared = curl_mom_eq_off+21 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\Theta\right)\right]_\theta\right)^2$ + Integer, Parameter :: curl_buoyancy_force_phi_squared = curl_mom_eq_off+22 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\Theta\right)\right]_\phi\right)^2$ + +! Integer, Parameter :: curl_buoyancy_pforce_theta = curl_mom_eq_off+23 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\overline{\Theta}\right)\right]_\theta$ +! Integer, Parameter :: curl_buoyancy_pforce_phi = curl_mom_eq_off+24 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\overline{\Theta}\right)\right]_\phi$ + +! Integer, Parameter :: curl_buoyancy_mforce_theta = curl_mom_eq_off+25 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\Theta'\right)\right]_\theta$ +! Integer, Parameter :: curl_buoyancy_mforce_phi = curl_mom_eq_off+26 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\Theta'\right)\right]_\phi$ - Integer, Parameter :: Coriolis_Force_r = curl_mom_eq_off+19 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\boldsymbol{v}\right]_r$ - Integer, Parameter :: Coriolis_Force_theta = curl_mom_eq_off+20 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\boldsymbol{v}\right]_\theta$ - Integer, Parameter :: Coriolis_Force_phi = curl_mom_eq_off+21 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\boldsymbol{v}\right]_\phi$ - Integer, Parameter :: Coriolis_pForce_r = curl_mom_eq_off+22 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\boldsymbol{v'}\right]_r$ - Integer, Parameter :: Coriolis_pForce_theta = curl_mom_eq_off+23 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\boldsymbol{v'}\right]_\theta$ - Integer, Parameter :: Coriolis_pForce_phi = curl_mom_eq_off+24 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\boldsymbol{v'}\right]_\phi$ + Integer, Parameter :: curl_coriolis_force_r = curl_mom_eq_off+27 ! :tex: $\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v})\right)\right]_r$ + Integer, Parameter :: curl_coriolis_force_theta = curl_mom_eq_off+28 ! :tex: $\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v})\right)\right]_\theta$ + Integer, Parameter :: curl_coriolis_force_phi = curl_mom_eq_off+29 ! :tex: $\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v})\right)\right]_\phi$ - Integer, Parameter :: Coriolis_mForce_r = curl_mom_eq_off+25 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\overline{\boldsymbol{v}}\right]_r$ - Integer, Parameter :: Coriolis_mForce_theta = curl_mom_eq_off+26 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\overline{\boldsymbol{v}}\right]_\theta$ - Integer, Parameter :: Coriolis_mForce_phi = curl_mom_eq_off+27 ! :tex: $-c_1\mathrm{f}_1\left[\boldsymbol{\hat{z}}\times\overline{\boldsymbol{v}}\right]_\phi$ + Integer, Parameter :: curl_coriolis_force_r_squared = curl_mom_eq_off+30 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v})\right)\right]_r\right)^2$ + Integer, Parameter :: curl_coriolis_force_theta_squared = curl_mom_eq_off+31 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v})\right)\right]_\theta\right)^2$ + Integer, Parameter :: curl_coriolis_force_phi_squared = curl_mom_eq_off+32 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v})\right)\right]_\phi\right)^2$ + +! Integer, Parameter :: curl_coriolis_pforce_r = curl_mom_eq_off+33 ! :tex: $\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v}')\right)\right]_r$ +! Integer, Parameter :: curl_coriolis_pforce_theta = curl_mom_eq_off+34 ! :tex: $\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v}')\right)\right]_\theta$ +! Integer, Parameter :: curl_coriolis_pforce_phi = curl_mom_eq_off+35 ! :tex: $\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v}')\right)\right]_\phi$ + +! Integer, Parameter :: curl_coriolis_mforce_r = curl_mom_eq_off+36 ! :tex: $\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\overline{\boldsymbol{\hat{v}}})\right)\right]_r$ +! Integer, Parameter :: curl_coriolis_mforce_theta = curl_mom_eq_off+37 ! :tex: $\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\overline{\boldsymbol{\hat{v}}})\right)\right]_\theta$ +! Integer, Parameter :: curl_coriolis_mforce_phi = curl_mom_eq_off+38 ! :tex: $\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\overline{\boldsymbol{\hat{v}}})\right)\right]_\phi$ ! Viscous forces - Integer, Parameter :: viscous_Force_r = curl_mom_eq_off+28 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D}}\right]_r$ - Integer, Parameter :: viscous_Force_theta = curl_mom_eq_off+29 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D}}\right]_\theta$ - Integer, Parameter :: viscous_Force_phi = curl_mom_eq_off+30 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D}}\right]_\phi$ + Integer, Parameter :: curl_viscous_force_r = curl_mom_eq_off+39 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_r$ + Integer, Parameter :: curl_viscous_force_theta = curl_mom_eq_off+40 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_\theta$ + Integer, Parameter :: curl_viscous_force_phi = curl_mom_eq_off+41 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_\phi$ - Integer, Parameter :: viscous_pForce_r = curl_mom_eq_off+31 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D'}}\right]_r$ - Integer, Parameter :: viscous_pForce_theta = curl_mom_eq_off+32 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D'}}\right]_\theta$ - Integer, Parameter :: viscous_pForce_phi = curl_mom_eq_off+33 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D'}}\right]_\phi$ +! Integer, Parameter :: curl_viscous_pforce_r = curl_mom_eq_off+42 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_r$ +! Integer, Parameter :: curl_viscous_pforce_theta = curl_mom_eq_off+43 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_\theta$ +! Integer, Parameter :: curl_viscous_pforce_phi = curl_mom_eq_off+44 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_\phi$ - Integer, Parameter :: viscous_mForce_r = curl_mom_eq_off+34 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal{D}}}\right]_r$ - Integer, Parameter :: viscous_mForce_theta = curl_mom_eq_off+35 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal{D}}}\right]_\theta$ - Integer, Parameter :: viscous_mForce_phi = curl_mom_eq_off+36 ! :tex: $c_5\left[\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal{D}}}\right]_\phi$ +! Integer, Parameter :: curl_viscous_mforce_r = curl_mom_eq_off+45 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_r$ +! Integer, Parameter :: curl_viscous_mforce_theta = curl_mom_eq_off+46 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\theta$ +! Integer, Parameter :: curl_viscous_mforce_phi = curl_mom_eq_off+47 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\phi$ ! Pressure forces - Integer, Parameter :: pressure_Force_r = curl_mom_eq_off+37 ! :tex: $ -c_3\mathrm{f}_1\frac{\partial}{\partial r}\left(\frac{P}{\mathrm{f}_1} \right)$ - Integer, Parameter :: pressure_Force_theta = curl_mom_eq_off+38 ! :tex: $ -c_3\frac{1}{r}\frac{\partial P}{\partial \theta}$ - Integer, Parameter :: pressure_Force_phi = curl_mom_eq_off+39 ! :tex: $ -c_3\frac{1}{r\mathrm{sin}\theta}\frac{\partial P}{\partial \phi}$ - - Integer, Parameter :: pressure_pForce_r = curl_mom_eq_off+40 ! :tex: $ -c_3\mathrm{f}_1\frac{\partial}{\partial r}\left(\frac{P'}{\mathrm{f}_1} \right)$ - Integer, Parameter :: pressure_pForce_theta = curl_mom_eq_off+41 ! :tex: $ -c_3\frac{1}{r}\frac{\partial P'}{\partial \theta}$ - Integer, Parameter :: pressure_pForce_phi = curl_mom_eq_off+42 ! :tex: $ -c_3\frac{1}{r\mathrm{sin}\theta}\frac{\partial P'}{\partial \phi}$ - - Integer, Parameter :: pressure_mForce_r = curl_mom_eq_off+43 ! :tex: $ -c_3\mathrm{f}_1\frac{\partial}{\partial r}\left(\frac{\overline{P}}{\mathrm{f}_1} \right)$ - Integer, Parameter :: pressure_mForce_theta = curl_mom_eq_off+44 ! :tex: $ -c_3\frac{1}{r}\frac{\partial \overline{P}}{\partial \theta}$ - Integer, Parameter :: pressure_mForce_phi = curl_mom_eq_off+45 ! :tex: $ -c_3\frac{1}{r\mathrm{sin}\theta}\frac{\partial \overline{P}}{\partial \phi}$ - - ! ell=0 pressure and buoyancy forces (r-direction only) - ! These are substracted out from the radial terms above - ! because they do not contribute to the dynamics. - ! The ell=0 buoyancy force is determined by the mean entropy profile. - ! The ell=0 pressure force is DEFINED such that it exactly balances the ell=0 buoyancy - Integer, Parameter :: buoyancy_force_ell0 = curl_mom_eq_off+46 ! :tex: $c_2\mathrm{f}_2\Theta_{00}$ - Integer, Parameter :: pressure_force_ell0_r = curl_mom_eq_off+47 ! :tex: $-c_3\mathrm{f}_1\frac{\partial}{\partial r}\left(\frac{P_{00}}{\mathrm{f}_1} \right)$ - - !/////////////////////////// Lorentz Forces /////////////////////////////// - ! ref%Lorentz_Coeff * (del x B) x B + Integer, Parameter :: curl_pressure_force_r = curl_mom_eq_off+48 ! :tex: + Integer, Parameter :: curl_pressure_force_theta = curl_mom_eq_off+49 ! :tex: + Integer, Parameter :: curl_pressure_force_phi = curl_mom_eq_off+50 ! :tex: + + Integer, Parameter :: curl_pressure_force_r_squared = curl_mom_eq_off+51 ! :tex: + Integer, Parameter :: curl_pressure_force_theta_squared = curl_mom_eq_off+52 ! :tex: + Integer, Parameter :: curl_pressure_force_phi_squared = curl_mom_eq_off+53 ! :tex: + +! Integer, Parameter :: curl_pressure_pforce_r = curl_mom_eq_off+54 ! :tex: +! Integer, Parameter :: curl_pressure_pforce_theta = curl_mom_eq_off+55 ! :tex: +! Integer, Parameter :: curl_pressure_pforce_phi = curl_mom_eq_off+56 ! :tex: + +! Integer, Parameter :: curl_pressure_mforce_r = curl_mom_eq_off+57 ! :tex: +! Integer, Parameter :: curl_pressure_mforce_theta = curl_mom_eq_off+58 ! :tex: +! Integer, Parameter :: curl_pressure_mforce_phi = curl_mom_eq_off+59 ! :tex: + + + !/////////////////////////// Lorentz forces /////////////////////////////// + ! ref%Lorentz_Coeff * del (del x B) x B ! ref%Lorentz_Coeff = 1/4pi when dimensional, Pr/(Pr_m E) when nondimesional ! j (below) is shorthand for ref%Lorentz_Coeff*delxB (not quite the current density) - Integer, Parameter :: j_cross_b_r = curl_mom_eq_off+48 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B}\right)\times\boldsymbol{B}\right]_r$ - Integer, Parameter :: j_cross_b_theta = curl_mom_eq_off+49 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B}\right)\times\boldsymbol{B}\right]_\theta$ - Integer, Parameter :: j_cross_b_phi = curl_mom_eq_off+50 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B}\right)\times\boldsymbol{B}\right]_\phi$ + Integer, Parameter :: curl_j_cross_b_r = curl_mom_eq_off+60 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_r$ + Integer, Parameter :: curl_j_cross_b_theta = curl_mom_eq_off+61 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\theta$ + Integer, Parameter :: curl_j_cross_b_phi = curl_mom_eq_off+62 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\phi$ + + Integer, Parameter :: curl_j_cross_b_r_squared = curl_mom_eq_off+63 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_r\right)^2$ + Integer, Parameter :: curl_j_cross_b_theta_squared = curl_mom_eq_off+64 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\theta\right)^2$ + Integer, Parameter :: curl_j_cross_b_phi_squared = curl_mom_eq_off+65 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\phi\right)^2$ + - Integer, Parameter :: jp_cross_bm_r = curl_mom_eq_off+51 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B'}\right)\times\overline{\boldsymbol{B}}\right]_r$ - Integer, Parameter :: jp_cross_bm_theta = curl_mom_eq_off+52 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B'}\right)\times\overline{\boldsymbol{B}}\right]_\theta$ - Integer, Parameter :: jp_cross_bm_phi = curl_mom_eq_off+53 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B'}\right)\times\overline{\boldsymbol{B}}\right]_\phi$ +! Integer, Parameter :: curl_jp_cross_bm_r = curl_mom_eq_off+66 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_r$ +! Integer, Parameter :: curl_jp_cross_bm_theta = curl_mom_eq_off+67 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_\theta$ +! Integer, Parameter :: curl_jp_cross_bm_phi = curl_mom_eq_off+68 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_\phi$ - Integer, Parameter :: jm_cross_bp_r = curl_mom_eq_off+54 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\overline{\boldsymbol{B}}\right)\times\boldsymbol{B'}\right]_r$ - Integer, Parameter :: jm_cross_bp_theta = curl_mom_eq_off+55 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\overline{\boldsymbol{B}}\right)\times\boldsymbol{B'}\right]_\theta$ - Integer, Parameter :: jm_cross_bp_phi = curl_mom_eq_off+56 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\overline{\boldsymbol{B}}\right)\times\boldsymbol{B'}\right]_\phi$ +! Integer, Parameter :: curl_jm_cross_bp_r = curl_mom_eq_off+69 ! :tex: $\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_r$ +! Integer, Parameter :: curl_jm_cross_bp_theta = curl_mom_eq_off+70 ! :tex: $\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_\theta$ +! Integer, Parameter :: curl_jm_cross_bp_phi = curl_mom_eq_off+71 ! :tex: $\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_\phi$ - Integer, Parameter :: jm_cross_bm_r = curl_mom_eq_off+57 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\overline{\boldsymbol{B}}\right)\times\overline{\boldsymbol{B}}\right]_r$ - Integer, Parameter :: jm_cross_bm_theta = curl_mom_eq_off+58 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\overline{\boldsymbol{B}}\right)\times\overline{\boldsymbol{B}}\right]_\theta$ - Integer, Parameter :: jm_cross_bm_phi = curl_mom_eq_off+59 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\overline{\boldsymbol{B}}\right)\times\overline{\boldsymbol{B}}\right]_\phi$ +! Integer, Parameter :: curl_jm_cross_bm_r = curl_mom_eq_off+72 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_r$ +! Integer, Parameter :: curl_jm_cross_bm_theta = curl_mom_eq_off+73 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_\theta$ +! Integer, Parameter :: curl_jm_cross_bm_phi = curl_mom_eq_off+74 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_\phi$ - Integer, Parameter :: jp_cross_bp_r = curl_mom_eq_off+60 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B'}\right)\times\boldsymbol{B'}\right]_r$ - Integer, Parameter :: jp_cross_bp_theta = curl_mom_eq_off+61 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B'}\right)\times\boldsymbol{B'}\right]_\theta$ - Integer, Parameter :: jp_cross_bp_phi = curl_mom_eq_off+62 ! :tex: $c_4\left[\left(\boldsymbol{\nabla}\times\boldsymbol{B'}\right)\times\boldsymbol{B'}\right]_\phi$ +! Integer, Parameter :: curl_jp_cross_bp_r = curl_mom_eq_off+74 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_r$ +! Integer, Parameter :: curl_jp_cross_bp_theta = curl_mom_eq_off+75 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_\theta$ +! Integer, Parameter :: curl_jp_cross_bp_phi = curl_mom_eq_off+76 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_\phi$ From a926426a4e6b60e27188c203cb74011b54311b13 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Fri, 17 Jul 2026 13:40:29 -0600 Subject: [PATCH 05/13] Refactor file (untested). --- src/Diagnostics/Diagnostics_Curl_Momentum.F90 | 91 +++++++++++++++---- 1 file changed, 73 insertions(+), 18 deletions(-) diff --git a/src/Diagnostics/Diagnostics_Curl_Momentum.F90 b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 index 317e6eb9e..af00d7632 100644 --- a/src/Diagnostics/Diagnostics_Curl_Momentum.F90 +++ b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 @@ -1,21 +1,48 @@ -Module Diagnostics_Custom +! +! Copyright (C) 2018 by the authors of the RAYLEIGH code. +! +! This file is part of RAYLEIGH. +! +! RAYLEIGH is free software; you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation; either version 3, or (at your option) +! any later version. +! +! RAYLEIGH is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with RAYLEIGH; see the file LICENSE. If not see +! . +! + +Module Diagnostics_Curl_Momentum Use Diagnostics_Base Implicit None + Contains - Subroutine Custom_MHD_Diagnostics(buffer) + Subroutine Compute_Curl_Momentum_Forces(buffer) + Real*8, Intent(InOut) :: buffer(1:,my_r%min:,my_theta%min:,1:) + Call Compute_Curl_Advection_Force(buffer) + Call Compute_Curl_Buoyancy_Force(buffer) + Call Compute_Curl_Magnetic_Force(buffer) + Call Compute_Curl_Coriolis_Force(buffer) + Call Compute_Curl_Pressure_Force(buffer) + End Subroutine + + Subroutine Compute_Curl_Advection_Force(buffer) Implicit None Real*8, Intent(InOut) :: buffer(1:,my_r%min:,my_theta%min:,1:) + Integer :: r, k, t + Real*8 :: pfactor(my_r%min:my_r%max) - Integer :: r,k, t pfactor(my_r%min:my_r%max) = ref%dpdr_w_term(my_r%min:my_r%max) & /ref%density(my_r%min:my_r%max) - - - !!!!!!!!!!!!!!!!!!!!!!!!! Advection Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - + !!!!!!!!!!!!!!!!!!!!!!!!! Advection Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! If (compute_quantity(curl_v_grad_v_r) .or. compute_quantity(curl_v_grad_v_r_squared)) Then If (compute_quantity(curl_v_grad_v_r)) Then @@ -53,7 +80,6 @@ Subroutine Custom_MHD_Diagnostics(buffer) Endif - If (compute_quantity(curl_v_grad_v_theta) .or. compute_quantity(curl_v_grad_v_theta_squared)) Then If (compute_quantity(curl_v_grad_v_theta)) Then DO_PSI @@ -127,8 +153,15 @@ Subroutine Custom_MHD_Diagnostics(buffer) Endif Endif + + End Subroutine Compute_Curl_Advection_Force + + Subroutine Compute_Curl_Buoyancy_Force(buffer) + Implicit None + Real*8, Intent(InOut) :: buffer(1:,my_r%min:,my_theta%min:,1:) + Integer :: r, k, t - !!!!!!!!!!!!!!!!!!!!!!!!!! Buoyancy Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!! Buoyancy Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! If (compute_quantity(curl_buoyancy_force_theta) .or. compute_quantity(curl_buoyancy_force_theta_squared)) Then @@ -167,8 +200,15 @@ Subroutine Custom_MHD_Diagnostics(buffer) Endif Endif - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Magnetic Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + End Subroutine Compute_Curl_Buoyancy_Force + + Subroutine Compute_Curl_Magnetic_Force(buffer) + Implicit None + Real*8, Intent(InOut) :: buffer(1:,my_r%min:,my_theta%min:,1:) + Integer :: r, k, t + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Magnetic Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! If (compute_quantity(curl_j_cross_b_r) .or. compute_quantity(curl_j_cross_b_r_squared)) Then If (compute_quantity(curl_j_cross_b_r)) Then @@ -286,7 +326,6 @@ Subroutine Custom_MHD_Diagnostics(buffer) Endif - If (compute_quantity(curl_j_cross_b_phi) .or. compute_quantity(curl_j_cross_b_phi_squared)) Then If (compute_quantity(curl_j_cross_b_phi)) Then DO_PSI @@ -340,7 +379,14 @@ Subroutine Custom_MHD_Diagnostics(buffer) Endif - !!!!!!!!!!!!!!!!!!!!!!!!!!!! Coriolis Force !!!!!!!!!!!!!!!!!!!!!!!!!! + End Subroutine Compute_Curl_Magnetic_Force + + Subroutine Compute_Curl_Coriolis_Force(buffer) + Implicit None + Real*8, Intent(InOut) :: buffer(1:,my_r%min:,my_theta%min:,1:) + Integer :: r, k, t + + !!!!!!!!!!!!!!!!!!!!!!!!!!!! Coriolis Force !!!!!!!!!!!!!!!!!!!!!!!!!! If (compute_quantity(curl_coriolis_force_r) .or. compute_quantity(curl_coriolis_force_r_squared)) Then If (compute_quantity(curl_coriolis_force_r)) Then @@ -404,12 +450,19 @@ Subroutine Custom_MHD_Diagnostics(buffer) Endif + End Subroutine Compute_Curl_Coriolis_Force + Subroutine Compute_Curl_Pressure_Force(buffer) + Implicit None + Real*8, Intent(InOut) :: buffer(1:,my_r%min:,my_theta%min:,1:) + Integer :: r, k, t + ! NOTE: pfactor is assumed constant in the derivation below + Real*8 :: pfactor(my_r%min:my_r%max) + pfactor(my_r%min:my_r%max) = ref%dpdr_w_term(my_r%min:my_r%max) & + /ref%density(my_r%min:my_r%max) !!!!!!!!!!!!!!!!!!!!!!!!!!! Pressure Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - If (compute_quantity(curl_pressure_force_r) .or. compute_quantity(curl_pressure_force_r_squared)) Then If (compute_quantity(curl_pressure_force_r)) Then DO_PSI @@ -467,5 +520,7 @@ Subroutine Custom_MHD_Diagnostics(buffer) Call Add_Quantity(qty) Endif Endif - - \ No newline at end of file + + End Subroutine Compute_Curl_Pressure_Force + +End Module Diagnostics_Curl_Momentum \ No newline at end of file From e974fbed2625fb1c29a34fc5cadb3c256eb5770b Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Sun, 19 Jul 2026 02:39:53 -0400 Subject: [PATCH 06/13] Implementation of the derivatives of the viscous forces. Untested but compiles. --- src/Diagnostics/Diagnostics_Base.F90 | 16 +- src/Diagnostics/Diagnostics_Curl_Momentum.F90 | 339 ++++++++++++++++-- src/Diagnostics/Diagnostics_Interface.F90 | 7 + src/Diagnostics/Diagnostics_Linear_Forces.F90 | 34 +- .../Diagnostics_Second_Derivatives.F90 | 1 + .../curl_momentum_equation_codes.F | 12 +- src/Include/indices.F | 1 + src/Makefile.fdeps | 6 +- src/object_list | 1 + 9 files changed, 367 insertions(+), 50 deletions(-) diff --git a/src/Diagnostics/Diagnostics_Base.F90 b/src/Diagnostics/Diagnostics_Base.F90 index 936842855..87cbbba97 100644 --- a/src/Diagnostics/Diagnostics_Base.F90 +++ b/src/Diagnostics/Diagnostics_Base.F90 @@ -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" @@ -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) @@ -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 diff --git a/src/Diagnostics/Diagnostics_Curl_Momentum.F90 b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 index af00d7632..e3d3d5681 100644 --- a/src/Diagnostics/Diagnostics_Curl_Momentum.F90 +++ b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 @@ -18,10 +18,17 @@ ! . ! +#include "indices.F" + Module Diagnostics_Curl_Momentum Use Diagnostics_Base + Use Spectral_Derivatives + Use Finite_Difference, Only : d_by_dx3d3 Implicit None + Integer, Allocatable :: vfdindmap(:,:) + Integer :: nvffields + Contains Subroutine Compute_Curl_Momentum_Forces(buffer) @@ -47,14 +54,18 @@ Subroutine Compute_Curl_Advection_Force(buffer) If (compute_quantity(curl_v_grad_v_r) .or. compute_quantity(curl_v_grad_v_r_squared)) Then If (compute_quantity(curl_v_grad_v_r)) Then DO_PSI - qty(PSI) = one_over_r(r) * ref%density(r) * (buffer(PSI,vr) * one_over_r(r) * DDBUFF(PSI,dvpdrdt) + & - buffer(PSI,dvrdt) * one_over_r(r) * buffer(PSI,dvpdr) + buffer(PSI,vtheta) * & - one_over_r(r) * DDBUFF(PSI,dvpdtdt) + buffer(PSI,dvtdt) * one_over_r(r) * & - buffer(PSI,dvpdt) + buffer(PSI,vphi) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvpdtdp) - & - csctheta(t) * csctheta(t) * costheta(t) * buffer(PSI,dvpdp) + buffer(PSI,dvrdt) + & - cottheta(t) * buffer(PSI,dvtdt) - csctheta(t) * & - csctheta(t) * buffer(PSI,vtheta)) + buffer(PSI,dvpdt) * one_over_r(r) * (csctheta(t) * & - buffer(PSI,dvpdp) + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) + & + qty(PSI) = one_over_r(r) * ref%density(r) * ( & + buffer(PSI,vr) * one_over_r(r) * DDBUFF(PSI,dvpdrdt) + & + buffer(PSI,dvrdt) * one_over_r(r) * buffer(PSI,dvpdr) + & + buffer(PSI,vtheta) * one_over_r(r) * DDBUFF(PSI,dvpdtdt) + & + buffer(PSI,dvtdt) * one_over_r(r) * buffer(PSI,dvpdt) + & + buffer(PSI,vphi) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvpdtdp) - & + csctheta(t) * csctheta(t) * costheta(t) * buffer(PSI,dvpdp) + & + buffer(PSI,dvrdt) + & + cottheta(t) * buffer(PSI,dvtdt) - & + csctheta(t) * csctheta(t) * buffer(PSI,vtheta)) + & + buffer(PSI,dvpdt) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvpdp) + & + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) + & one_over_r(r) * costheta(t) * csctheta(t) * ref%density(r) * (buffer(PSI,vr) * & one_over_r(r) * buffer(PSI,dvpdr) + buffer(PSI,vtheta) * one_over_r(r) * & buffer(PSI,dvpdt) + buffer(PSI,vphi) * one_over_r(r) * (csctheta(t) * & @@ -92,16 +103,20 @@ Subroutine Compute_Curl_Advection_Force(buffer) one_over_r(r) * ref%density(r) * (buffer(PSI,vr) * one_over_r(r) * buffer(PSI,dvpdr) + & buffer(PSI,vtheta) * one_over_r(r) * buffer(PSI,dvpdt) + buffer(PSI,vphi) * one_over_r(r) * & (csctheta(t) * buffer(PSI,dvpdp) + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) - & - ref%density(r) * (buffer(PSI,vr) * one_over_r(r) * DDBUFF(PSI,dvpdrdr) + buffer(PSI,dvrdr) * & - one_over_r(r) * buffer(PSI,dvpdr) - buffer(PSI,vr) * one_over_r(r) * one_over_r(r) * & - buffer(PSI,dvpdr) + buffer(PSI,vtheta) * one_over_r(r) * DDBUFF(PSI,dvpdrdt) + & - buffer(PSI,dvtdr) * one_over_r(r) * buffer(PSI,dvpdt) - buffer(PSI,vtheta) * one_over_r(r) * & - one_over_r(r) * buffer(PSI,dvpdt) + & - buffer(PSI,vtheta) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvpdrdp) + & - buffer(PSI,dvrdr) + cottheta(t) * buffer(PSI,dvtdr)) + & - buffer(PSI,dvpdr) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvpdp) + buffer(PSI,vr) + & - cottheta(t) * buffer(PSI,vtheta)) - buffer(PSI,vphi) * one_over_r(r) * one_over_r(r) * & - (csctheta(t) * buffer(PSI,dvpdp) + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) - & + ref%density(r) * (buffer(PSI,vr) * one_over_r(r) * DDBUFF(PSI,dvpdrdr) + & + buffer(PSI,dvrdr) * one_over_r(r) * buffer(PSI,dvpdr) - & + buffer(PSI,vr) * one_over_r(r) * one_over_r(r) * buffer(PSI,dvpdr) + & + buffer(PSI,vtheta) * one_over_r(r) * DDBUFF(PSI,dvpdrdt) + & + buffer(PSI,dvtdr) * one_over_r(r) * buffer(PSI,dvpdt) - & + buffer(PSI,vtheta) * one_over_r(r) * one_over_r(r) * buffer(PSI,dvpdt) + & + buffer(PSI,vtheta) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvpdrdp) + & + buffer(PSI,dvrdr) + & + cottheta(t) * buffer(PSI,dvtdr)) + & + buffer(PSI,dvpdr) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvpdp) + & + buffer(PSI,vr) + & + cottheta(t) * buffer(PSI,vtheta)) - & + buffer(PSI,vphi) * one_over_r(r) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvpdp) + & + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) - & ref%dlnrho(r) * & ref%density(r) * (buffer(PSI,vr) * one_over_r(r) * buffer(PSI,dvpdr) + & buffer(PSI,vtheta) * one_over_r(r) * buffer(PSI,dvpdt) + buffer(PSI,vphi) * one_over_r(r) * & @@ -121,25 +136,34 @@ Subroutine Compute_Curl_Advection_Force(buffer) If (compute_quantity(curl_v_grad_v_phi) .or. compute_quantity(curl_v_grad_v_phi_squared)) Then If (compute_quantity(curl_v_grad_v_theta)) Then DO_PSI - qty(PSI) = ref%density(r) * (buffer(PSI,vr) * DDBUFF(PSI,dvtdrdr) + buffer(PSI,dvrdr) * buffer(PSI,dvtdr) + & - buffer(PSI,vtheta) * one_over_r(r) * (DDBUFF(PSI,dvtdrdt) + buffer(PSI,dvrdr)) + & - buffer(PSI,dvtdr) * one_over_r(r) * (buffer(PSI,dvtdt) + buffer(PSI,vr)) - & - buffer(PSI,vtheta) * one_over_r(r) * one_over_r(r) * (buffer(PSI,dvtdt) + buffer(PSI,vr)) + & - buffer(PSI,vtheta) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvtdrdp) - cottheta(t) * & - buffer(PSI,dvpdr)) + buffer(PSI,dvpdr) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvtdp) - & - cottheta(t) * buffer(PSI,vphi)) - buffer(PSI,vphi) * one_over_r(r) * one_over_r(r) * & - (csctheta(t) * buffer(PSI,dvtdp) - cottheta(t) * buffer(PSI,vphi))) + & + qty(PSI) = ref%density(r) * (buffer(PSI,vr) * DDBUFF(PSI,dvtdrdr) + & + buffer(PSI,dvrdr) * buffer(PSI,dvtdr) + & + buffer(PSI,vtheta) * one_over_r(r) * (DDBUFF(PSI,dvtdrdt) + & + buffer(PSI,dvrdr)) + & + buffer(PSI,dvtdr) * one_over_r(r) * (buffer(PSI,dvtdt) + & + buffer(PSI,vr)) - & + buffer(PSI,vtheta) * one_over_r(r) * one_over_r(r) * (buffer(PSI,dvtdt) + & + buffer(PSI,vr)) + & + buffer(PSI,vtheta) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvtdrdp) - & + cottheta(t) * buffer(PSI,dvpdr)) + & + buffer(PSI,dvpdr) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvtdp) - & + cottheta(t) * buffer(PSI,vphi)) - & + buffer(PSI,vphi) * one_over_r(r) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvtdp) - & + cottheta(t) * buffer(PSI,vphi))) + & ref%dlnrho(r) * ref%density(r) * (buffer(PSI,vr) * buffer(PSI,dvtdr) + buffer(PSI,vtheta) * & one_over_r(r) * (buffer(PSI,dvtdt) + buffer(PSI,vr)) + buffer(PSI,vphi) * one_over_r(r) * & (csctheta(t) * buffer(PSI,dvtdp) - cottheta(t) * buffer(PSI,vphi))) + & one_over_r(r) * ref%density(r) * (buffer(PSI,vr) * buffer(PSI,dvtdr) + buffer(PSI,vtheta) * & one_over_r(r) * (buffer(PSI,dvtdt) + buffer(PSI,vr)) + buffer(PSI,vphi) * one_over_r(r) * & (csctheta(t) * buffer(PSI,dvtdp) - cottheta(t) * buffer(PSI,vphi))) - & - one_over_r(r) * ref%density(r) * (buffer(PSI,vr) * DDBUFF(PSI,dvrdrdt) + buffer(PSI,dvrdt) * & - buffer(PSI,dvrdr) + buffer(PSI,vtheta) * one_over_r(r) * (DDBUFF(PSI,dvrdtdt) + buffer(PSI,dvtdt)) + & + one_over_r(r) * ref%density(r) * (buffer(PSI,vr) * DDBUFF(PSI,dvrdrdt) + & + buffer(PSI,dvrdt) * buffer(PSI,dvrdr) + buffer(PSI,vtheta) * one_over_r(r) * & + (DDBUFF(PSI,dvrdtdt) + buffer(PSI,dvtdt)) + & buffer(PSI,dvtdt) * one_over_r(r) * (buffer(PSI,dvrdt) - buffer(PSI,vtheta)) + & - buffer(PSI,vphi) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvrdtdp) - csctheta(t) * csctheta(t) * & - costheta(t) * buffer(PSI,dvrdp) - buffer(PSI,dvpdt)) + buffer(PSI,dvpdt) * one_over_r(r) * & + buffer(PSI,vphi) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvrdtdp) - & + csctheta(t) * csctheta(t) * costheta(t) * buffer(PSI,dvrdp) - & + buffer(PSI,dvpdt)) + & + buffer(PSI,dvpdt) * one_over_r(r) * & (csctheta(t) * buffer(PSI,dvrdp) - buffer(PSI,vphi))) END_DO @@ -221,8 +245,10 @@ Subroutine Compute_Curl_Magnetic_Force(buffer) buffer(PSI,bphi) + & one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,btheta) * & buffer(PSI,dbpdt) - & - one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * buffer(PSI,dbtdt) * buffer(PSI,dbtdp) - & - one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * buffer(PSI,btheta) * DDBUFF(PSI,dbtdtdp) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * & + buffer(PSI,dbtdt) * buffer(PSI,dbtdp) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * & + buffer(PSI,btheta) * DDBUFF(PSI,dbtdtdp) + & one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * buffer(PSI,btheta) * & buffer(PSI,dbtdp) + & one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * & @@ -522,5 +548,252 @@ Subroutine Compute_Curl_Pressure_Force(buffer) Endif End Subroutine Compute_Curl_Pressure_Force + + Subroutine Initialize_Grad_Viscous_Force() + Implicit None + integer :: nvfind, nvfdind, vfoff + integer :: nvfdrfields, nvfdtfields, nvfdpfields + integer :: vfdfcount(3,2) ! buffer sizes + Integer :: vf_i(9) ! indices to vforce_buffer + Logical :: compute_vforce_i_dj(9,3) + Integer :: i, j + + dvf_r_dt = -1 + dvf_r_dp = -1 + dvf_t_dr = -1 + dvf_t_dp = -1 + dvf_p_dr = -1 + dvf_p_dt = -1 + dvfp_r_dt = -1 + dvfp_r_dp = -1 + dvfp_t_dr = -1 + dvfp_t_dp = -1 + dvfp_p_dr = -1 + dvfp_p_dt = -1 + dvfm_r_dt = -1 + dvfm_r_dp = -1 + dvfm_t_dr = -1 + dvfm_t_dp = -1 + dvfm_p_dr = -1 + dvfm_p_dt = -1 + + vf_i = [vf_r, vf_t, vf_p, vfp_r, vfp_t, vfp_p, vfm_r, vfm_t, vfm_p] + compute_vforce_i_dj = .false. + + vfoff = 0 + If (sometimes_compute(curl_viscous_force_r)) Then + compute_vforce_i_dj(2,3) = .true. + compute_vforce_i_dj(3,2) = .true. + Endif + + If (sometimes_compute(curl_viscous_force_theta)) Then + compute_vforce_i_dj(1,3) = .true. + compute_vforce_i_dj(3,1) = .true. + Endif + + If (sometimes_compute(curl_viscous_force_phi)) Then + compute_vforce_i_dj(1,2) = .true. + compute_vforce_i_dj(2,1) = .true. + Endif + + vfoff = 3 + If (sometimes_compute(curl_viscous_pforce_r)) Then + compute_vforce_i_dj(vfoff+2,3) = .true. + compute_vforce_i_dj(vfoff+3,2) = .true. + Endif + + If (sometimes_compute(curl_viscous_pforce_theta)) Then + compute_vforce_i_dj(vfoff+1,3) = .true. + compute_vforce_i_dj(vfoff+3,1) = .true. + Endif + + If (sometimes_compute(curl_viscous_pforce_phi)) Then + compute_vforce_i_dj(vfoff+1,2) = .true. + compute_vforce_i_dj(vfoff+2,1) = .true. + Endif + + vfoff = 6 + If (sometimes_compute(curl_viscous_mforce_r)) Then + compute_vforce_i_dj(vfoff+2,3) = .true. + compute_vforce_i_dj(vfoff+3,2) = .true. + Endif + + If (sometimes_compute(curl_viscous_mforce_theta)) Then + compute_vforce_i_dj(vfoff+1,3) = .true. + compute_vforce_i_dj(vfoff+3,1) = .true. + Endif + + If (sometimes_compute(curl_viscous_mforce_phi)) Then + compute_vforce_i_dj(vfoff+1,2) = .true. + compute_vforce_i_dj(vfoff+2,1) = .true. + Endif + + ! work out how many vf fields we'll be taking the derivative of + nvffields = count(count(compute_vforce_i_dj, dim=2) .gt. 0) + Allocate(vfdindmap(nvffields,4)) + vfdindmap(:,:) = -1 + need_vforce_derivatives = nvffields .gt. 0 + + ! next assign indices to vf_i and vforce derivative entries + ! this loop is designed to make sure the new derivative fields are indexed + ! in r, t, p order so that we can grow the buffers appropriately + nvfdind = 0 + do j = 1, 3 + nvfind = 0 + do i = 1, 9 + if (count(compute_vforce_i_dj(i,:)) .gt. 0) then + nvfind = nvfind + 1 + ! assign indices to vforce_buffer in first column (on first outer loop) + if (j .eq. 1) vfdindmap(nvfind, 1) = vf_i(i) + if (compute_vforce_i_dj(i,j)) then + nvfdind = nvfdind + 1 + if ((i .eq. 1) .and. (j .eq. 2)) dvf_r_dt = nvfdind + if ((i .eq. 1) .and. (j .eq. 3)) dvf_r_dp = nvfdind + if ((i .eq. 2) .and. (j .eq. 1)) dvf_t_dr = nvfdind + if ((i .eq. 2) .and. (j .eq. 3)) dvf_t_dp = nvfdind + if ((i .eq. 3) .and. (j .eq. 1)) dvf_p_dr = nvfdind + if ((i .eq. 3) .and. (j .eq. 2)) dvf_p_dt = nvfdind + vfoff = 3 + if ((i .eq. vfoff+1) .and. (j .eq. 2)) dvfp_r_dt = nvfdind + if ((i .eq. vfoff+1) .and. (j .eq. 3)) dvfp_r_dp = nvfdind + if ((i .eq. vfoff+2) .and. (j .eq. 1)) dvfp_t_dr = nvfdind + if ((i .eq. vfoff+2) .and. (j .eq. 3)) dvfp_t_dp = nvfdind + if ((i .eq. vfoff+3) .and. (j .eq. 1)) dvfp_p_dr = nvfdind + if ((i .eq. vfoff+3) .and. (j .eq. 2)) dvfp_p_dt = nvfdind + vfoff = 6 + if ((i .eq. vfoff+1) .and. (j .eq. 2)) dvfm_r_dt = nvfdind + if ((i .eq. vfoff+1) .and. (j .eq. 3)) dvfm_r_dp = nvfdind + if ((i .eq. vfoff+2) .and. (j .eq. 1)) dvfm_t_dr = nvfdind + if ((i .eq. vfoff+2) .and. (j .eq. 3)) dvfm_t_dp = nvfdind + if ((i .eq. vfoff+3) .and. (j .eq. 1)) dvfm_p_dr = nvfdind + if ((i .eq. vfoff+3) .and. (j .eq. 2)) dvfm_p_dt = nvfdind + vfdindmap(nvfind, j+1) = nvfdind + endif + endif + enddo + enddo + + ! work out how many of each type of derivative we're taking + nvfdrfields = count(compute_vforce_i_dj(:,1)) + nvfdtfields = count(compute_vforce_i_dj(:,2)) + nvfdpfields = count(compute_vforce_i_dj(:,3)) + + ! size the buffers at each config stage + vfdfcount(1,1) = nvffields + nvfdrfields ! config 1a + vfdfcount(2,1) = nvffields + nvfdrfields + nvfdtfields ! 2a + vfdfcount(3,1) = nvffields + nvfdrfields + nvfdtfields + nvfdpfields ! 3a + vfdfcount(3,2) = nvffields ! 3b + vfdfcount(2,2) = nvffields ! 2b + vfdfcount(1,2) = nvffields + nvfdrfields ! 1b + + Call d_vforce_buffer%init(field_count = vfdfcount, config = 'p3b') + + Call d_vforce_buffer%construct('p3a') + + Call d_vforce_buffer%deconstruct('p3a') + + End Subroutine Initialize_Grad_Viscous_Force + + Subroutine Grad_Viscous_Force() + Implicit None + Integer :: i, r, k, t + + call d_vforce_buffer%construct('p3b') + d_vforce_buffer%config = 'p3b' + + ! load the fields we want to take derivatives of + do i = 1, nvffields + d_vforce_buffer%p3b(:,:,:,i) = vforce_buffer(:,:,:,vfdindmap(i, 1)) + enddo + + ! transform to Fourier m space + Call fft_to_spectral(d_vforce_buffer%p3b, rsc = .true.) + + ! reform to hybrid rlm space + call d_vforce_buffer%reform() ! move to p2b + + ! deallocate p3b + call d_vforce_buffer%deconstruct('p3b') + + ! allocate spectral buffer and transform + call d_vforce_buffer%construct('s2b') + call Legendre_Transform(d_vforce_buffer%p2b, d_vforce_buffer%s2b) + + ! deallocate p2b + call d_vforce_buffer%deconstruct('p2b') + + ! reform + call d_vforce_buffer%reform() ! move to p1b + + ! do a little gymnastics with p1a and p1b + call d_vforce_buffer%construct('p1a') + if (chebyshev) then + ! store chebyshev coefficients in p1a and dealias + call gridcp%to_Spectral(d_vforce_buffer%p1b, d_vforce_buffer%p1a) + call gridcp%dealias_buffer(d_vforce_buffer%p1a) + else + d_vforce_buffer%p1a = d_vforce_buffer%p1b + end if + + ! take d_by_dr + ! (and transform back to grid space if in Chebyshev) + if (chebyshev) then + do i = 1, nvffields + if (vfdindmap(i,2) .gt. 0) then + call gridcp%d_by_dr_cp(i, vfdindmap(i,2), d_vforce_buffer%p1a, 1) + endif + enddo + call gridcp%from_spectral(d_vforce_buffer%p1a, d_vforce_buffer%p1b) + ! FIXME: different to second derivatives + d_vforce_buffer%p1a = d_vforce_buffer%p1b + else + do i = 1, nvffields + if (vfdindmap(i,2) .gt. 0) then + call d_by_dx3d3(i, vfdindmap(i,2), d_vforce_buffer%p1a,1) + endif + enddo + end if + + ! moving back, but first we need grid space in p1a + d_vforce_buffer%config = 'p1a' + call d_vforce_buffer%deconstruct('p1b') + + ! reform and start moving back + call d_vforce_buffer%reform() ! now in s2a + + ! take theta derivatives + do i = 1, nvffields + if (vfdindmap(i,3) .gt. 0) then + call d_by_dtheta(d_vforce_buffer%s2a, i, vfdindmap(i,3)) + endif + enddo + + call d_vforce_buffer%construct('p2a') + call Legendre_Transform(d_vforce_buffer%s2a, d_vforce_buffer%p2a) + call d_vforce_buffer%deconstruct('s2a') + + ! reform + call d_vforce_buffer%reform() ! move to p2a + + ! take d_by_dphi derivatives + do i = 1, nvffields + if (vfdindmap(i,4) .gt. 0) then + call d_by_dphi(d_vforce_buffer%p3a, i, vfdindmap(i,4)) + endif + enddo + + ! transform to grid space + call FFT_To_Physical(d_vforce_buffer%p3a, rsc=.true.) + + ! Convert sintheta*{dxdt} to dxdt + do i = 1, nvffields + if (vfdindmap(i,3) .gt. 0) then + DO_PSI + d_vforce_buffer%p3a(PSI,i) = d_vforce_buffer%p3a(PSI,i)*csctheta(t) + END_DO + end if + enddo + + End Subroutine Grad_Viscous_Force -End Module Diagnostics_Curl_Momentum \ No newline at end of file +End Module Diagnostics_Curl_Momentum diff --git a/src/Diagnostics/Diagnostics_Interface.F90 b/src/Diagnostics/Diagnostics_Interface.F90 index 86469b137..1fe54e1c0 100755 --- a/src/Diagnostics/Diagnostics_Interface.F90 +++ b/src/Diagnostics/Diagnostics_Interface.F90 @@ -31,6 +31,7 @@ Module Diagnostics_Interface Use Diagnostics_Base Use Diagnostics_Second_Derivatives + Use Diagnostics_Curl_Momentum Use Diagnostics_Mean_Correction @@ -189,6 +190,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 @@ -338,6 +343,8 @@ Subroutine Initialize_Diagnostics() Call Initialize_Second_Derivatives() + Call Initialize_Grad_Viscous_Force() + Call Initialize_Diagnostics_Buffer() End Subroutine Initialize_Diagnostics diff --git a/src/Diagnostics/Diagnostics_Linear_Forces.F90 b/src/Diagnostics/Diagnostics_Linear_Forces.F90 index e4598baed..e93a5348d 100644 --- a/src/Diagnostics/Diagnostics_Linear_Forces.F90 +++ b/src/Diagnostics/Diagnostics_Linear_Forces.F90 @@ -519,56 +519,74 @@ Subroutine Initialize_Viscous_Force() If (compute_quantity(viscous_force_r) .or. & compute_quantity(visc_work) .or. & + compute_quantity(curl_viscous_force_theta) .or. & + compute_quantity(curl_viscous_force_phi) .or. & compute_quantity(viscous_mforce_r)) Then nvf = nvf+1 vf_r = nvf Endif 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_phi)) Then nvf = nvf + 1 vf_t = nvf Endif 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_theta)) Then nvf = nvf + 1 vf_p = nvf Endif 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 nvf = nvf+1 vfp_r = nvf Endif 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 nvf = nvf + 1 vfp_t = nvf Endif 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 nvf = nvf + 1 vfp_p = nvf Endif 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 nvf = nvf+1 vfm_r = nvf Endif 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 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 + compute_quantity(visc_work_mm).or. & + compute_quantity(curl_viscous_mforce_r) .or. & + compute_quantity(curl_viscous_mforce_theta)) Then nvf = nvf + 1 vfm_p = nvf Endif diff --git a/src/Diagnostics/Diagnostics_Second_Derivatives.F90 b/src/Diagnostics/Diagnostics_Second_Derivatives.F90 index 27ae6b8ec..bc752f0cf 100644 --- a/src/Diagnostics/Diagnostics_Second_Derivatives.F90 +++ b/src/Diagnostics/Diagnostics_Second_Derivatives.F90 @@ -400,6 +400,7 @@ Subroutine Compute_Second_Derivatives(inbuffer) Call d_by_dx3d3(i,j,d2buffer%p1a,1) Enddo Endif + ! CRW: Is this right in finite differences? d2buffer%p1a=d2buffer%p1b Call d2buffer%deconstruct('p1b') diff --git a/src/Diagnostics/curl_momentum_equation_codes.F b/src/Diagnostics/curl_momentum_equation_codes.F index b6180face..9cbf86669 100644 --- a/src/Diagnostics/curl_momentum_equation_codes.F +++ b/src/Diagnostics/curl_momentum_equation_codes.F @@ -85,13 +85,13 @@ Integer, Parameter :: curl_viscous_force_theta = curl_mom_eq_off+40 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_\theta$ Integer, Parameter :: curl_viscous_force_phi = curl_mom_eq_off+41 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_\phi$ -! Integer, Parameter :: curl_viscous_pforce_r = curl_mom_eq_off+42 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_r$ -! Integer, Parameter :: curl_viscous_pforce_theta = curl_mom_eq_off+43 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_\theta$ -! Integer, Parameter :: curl_viscous_pforce_phi = curl_mom_eq_off+44 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_\phi$ + Integer, Parameter :: curl_viscous_pforce_r = curl_mom_eq_off+42 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_r$ + Integer, Parameter :: curl_viscous_pforce_theta = curl_mom_eq_off+43 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_\theta$ + Integer, Parameter :: curl_viscous_pforce_phi = curl_mom_eq_off+44 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_\phi$ -! Integer, Parameter :: curl_viscous_mforce_r = curl_mom_eq_off+45 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_r$ -! Integer, Parameter :: curl_viscous_mforce_theta = curl_mom_eq_off+46 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\theta$ -! Integer, Parameter :: curl_viscous_mforce_phi = curl_mom_eq_off+47 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\phi$ + Integer, Parameter :: curl_viscous_mforce_r = curl_mom_eq_off+45 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_r$ + Integer, Parameter :: curl_viscous_mforce_theta = curl_mom_eq_off+46 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\theta$ + Integer, Parameter :: curl_viscous_mforce_phi = curl_mom_eq_off+47 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\phi$ ! Pressure forces Integer, Parameter :: curl_pressure_force_r = curl_mom_eq_off+48 ! :tex: diff --git a/src/Include/indices.F b/src/Include/indices.F index 0b0abcdb2..c5893e5fe 100644 --- a/src/Include/indices.F +++ b/src/Include/indices.F @@ -35,4 +35,5 @@ #define PSI k,r,t #define PSI2 r,t #define DDBUFF d2buffer%p3a +#define VFDBUFF d_vforce_buffer%p3a diff --git a/src/Makefile.fdeps b/src/Makefile.fdeps index f3c80d6fb..c67f7a1f5 100644 --- a/src/Makefile.fdeps +++ b/src/Makefile.fdeps @@ -10,13 +10,14 @@ Controls.o : Controls.F90 BufferedOutput.o Diagnostics_ADotGradB.o : Diagnostics_ADotGradB.F90 indices.F Diagnostics_Base.o Diagnostics_Angular_Momentum.o : Diagnostics_Angular_Momentum.F90 indices.F Diagnostics_Base.o Diagnostics_Axial_Field.o : Diagnostics_Axial_Field.F90 indices.F Diagnostics_Base.o -Diagnostics_Base.o : Diagnostics_Base.F90 scalars_field_codes.F axial_field_codes.F turbKE_codes.F me_equation_codes.F ke_equation_codes.F amom_equation_codes.F induction_equation_codes.F thermal_equation_codes.F momentum_equation_codes.F magnetic_energy_codes.F current_density_codes.F magnetic_field_codes.F thermal_energy_codes.F thermal_field_codes.F kinetic_energy_codes.F vorticity_field_codes.F mass_flux_codes.F velocity_field_codes.F indices.F PDE_Coefficients.o Math_Constants.o Fields.o Spherical_IO.o ProblemSize.o +Diagnostics_Base.o : Diagnostics_Base.F90 scalars_field_codes.F axial_field_codes.F turbKE_codes.F me_equation_codes.F ke_equation_codes.F amom_equation_codes.F induction_equation_codes.F thermal_equation_codes.F curl_momentum_equation_codes.F momentum_equation_codes.F magnetic_energy_codes.F current_density_codes.F magnetic_field_codes.F thermal_energy_codes.F thermal_field_codes.F kinetic_energy_codes.F vorticity_field_codes.F mass_flux_codes.F velocity_field_codes.F indices.F PDE_Coefficients.o Math_Constants.o Fields.o Spherical_IO.o ProblemSize.o +Diagnostics_Curl_Momentum.o : Diagnostics_Curl_Momentum.F90 indices.F Finite_Difference.o Spectral_Derivatives.o Diagnostics_Base.o Diagnostics_Current_Density.o : Diagnostics_Current_Density.F90 indices.F Diagnostics_Base.o Diagnostics_Custom.o : Diagnostics_Custom.F90 indices.F Diagnostics_Base.o Diagnostics_Energies.o : Diagnostics_Energies.F90 indices.F Diagnostics_Base.o Diagnostics_Induction.o : Diagnostics_Induction.F90 indices.F Diagnostics_ADotGradB.o Diagnostics_Base.o Diagnostics_Inertial_Forces.o : Diagnostics_Inertial_Forces.F90 indices.F Diagnostics_ADotGradB.o Diagnostics_Base.o -Diagnostics_Interface.o : Diagnostics_Interface.F90 indices.F Diagnostics_Scalars.o Diagnostics_Custom.o Diagnostics_Miscellaneous.o Diagnostics_Poynting_Flux.o Diagnostics_Induction.o Diagnostics_Axial_Field.o Diagnostics_TurbKE_Budget.o Diagnostics_KE_Flux.o Diagnostics_Lorentz_Forces.o Diagnostics_Angular_Momentum.o Diagnostics_Inertial_Forces.o Diagnostics_Linear_Forces.o Diagnostics_Current_Density.o Diagnostics_Vorticity_Field.o Diagnostics_Thermal_Equation.o Diagnostics_Thermal_Energies.o Diagnostics_Thermodynamic_Gradients.o Diagnostics_Energies.o Diagnostics_Magnetic_Field.o Diagnostics_Velocity_Field.o Diagnostics_Mean_Correction.o Diagnostics_Second_Derivatives.o Diagnostics_Base.o Math_Constants.o PDE_Coefficients.o Legendre_Polynomials.o Fields.o Spherical_IO.o Controls.o ProblemSize.o +Diagnostics_Interface.o : Diagnostics_Interface.F90 indices.F Diagnostics_Scalars.o Diagnostics_Custom.o Diagnostics_Miscellaneous.o Diagnostics_Poynting_Flux.o Diagnostics_Induction.o Diagnostics_Axial_Field.o Diagnostics_TurbKE_Budget.o Diagnostics_KE_Flux.o Diagnostics_Lorentz_Forces.o Diagnostics_Angular_Momentum.o Diagnostics_Inertial_Forces.o Diagnostics_Linear_Forces.o Diagnostics_Current_Density.o Diagnostics_Vorticity_Field.o Diagnostics_Thermal_Equation.o Diagnostics_Thermal_Energies.o Diagnostics_Thermodynamic_Gradients.o Diagnostics_Energies.o Diagnostics_Magnetic_Field.o Diagnostics_Velocity_Field.o Diagnostics_Mean_Correction.o Diagnostics_Curl_Momentum.o Diagnostics_Second_Derivatives.o Diagnostics_Base.o Math_Constants.o PDE_Coefficients.o Legendre_Polynomials.o Fields.o Spherical_IO.o Controls.o ProblemSize.o Diagnostics_KE_Flux.o : Diagnostics_KE_Flux.F90 indices.F Diagnostics_Base.o Diagnostics_Linear_Forces.o : Diagnostics_Linear_Forces.F90 indices.F Diagnostics_Base.o Diagnostics_Lorentz_Forces.o : Diagnostics_Lorentz_Forces.F90 indices.F Diagnostics_Base.o @@ -79,6 +80,7 @@ Timers.o : Timers.F90 SendReceive.o Parallel_Framework.o Timing.o Timing.o : Timing.F90 Ra_MPI_Base.o amom_equation_codes.o : amom_equation_codes.F axial_field_codes.o : axial_field_codes.F +curl_momentum_equation_codes.o : curl_momentum_equation_codes.F current_density_codes.o : current_density_codes.F indices.o : indices.F induction_equation_codes.o : induction_equation_codes.F diff --git a/src/object_list b/src/object_list index 4c52fdc14..7901e0f4d 100755 --- a/src/object_list +++ b/src/object_list @@ -16,6 +16,7 @@ POBJ = Controls.o ClockInfo.o ProblemSize.o PDE_Coefficients.o Fields.o \ Diagnostics_Current_Density.o Diagnostics_Vorticity_Field.o Diagnostics_Axial_Field.o \ Diagnostics_Thermal_Energies.o Diagnostics_Thermal_Equation.o\ Diagnostics_Thermodynamic_Gradients.o Diagnostics_Linear_Forces.o \ + Diagnostics_Curl_Momentum.o \ Diagnostics_Energies.o \ Diagnostics_Lorentz_Forces.o Diagnostics_Induction.o \ Diagnostics_Inertial_Forces.o Diagnostics_Angular_Momentum.o \ From ceaa2248a4eca3f921746b4abd50ee28105d397c Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Sun, 19 Jul 2026 16:11:06 +0000 Subject: [PATCH 07/13] Trying to update docs to reflect new curl of momentum quantity codes. This involved adding some dummy ones as blank tex is not allowed. The parser also doesn't seem to recognize commented out lines. --- doc/build_menu_pdf.py | 2 + .../curl_momentum_equation.rst | 82 +++++++++++++++++++ doc/source/diagnostic_codes/qcodes.rst | 1 + .../curl_momentum_equation_codes.F | 24 +++--- 4 files changed, 97 insertions(+), 12 deletions(-) create mode 100644 doc/source/diagnostic_codes/curl_momentum_equation.rst diff --git a/doc/build_menu_pdf.py b/doc/build_menu_pdf.py index 278711417..402b93a08 100644 --- a/doc/build_menu_pdf.py +++ b/doc/build_menu_pdf.py @@ -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]) @@ -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') diff --git a/doc/source/diagnostic_codes/curl_momentum_equation.rst b/doc/source/diagnostic_codes/curl_momentum_equation.rst new file mode 100644 index 000000000..0729d2762 --- /dev/null +++ b/doc/source/diagnostic_codes/curl_momentum_equation.rst @@ -0,0 +1,82 @@ +Curl Momentum Equation +==================================================================== + +===================================================================================================================================================================== ====== ======================================= + :math:`\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_r` 1301 curl\_v\_grad\_v\_r + :math:`\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\theta` 1302 curl\_v\_grad\_v\_theta + :math:`\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\phi` 1303 curl\_v\_grad\_v\_phi + :math:`\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_r\right)^2` 1304 curl\_v\_grad\_v\_r\_squared + :math:`\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\theta\right)^2` 1305 curl\_v\_grad\_v\_theta\_squared + :math:`\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\phi\right)^2` 1306 curl\_v\_grad\_v\_phi\_squared + :math:`\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_r` 1307 curl\_vp\_grad\_vm\_r + :math:`\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\theta` 1308 curl\_vp\_grad\_vm\_theta + :math:`\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\phi` 1309 curl\_vp\_grad\_vm\_phi + :math:`\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_r` 1310 curl\_vm\_grad\_vp\_r + :math:`\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\theta` 1311 curl\_vm\_grad\_vp\_theta + :math:`\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\phi` 1312 curl\_vm\_grad\_vp\_phi + :math:`\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_r` 1313 curl\_vp\_grad\_vp\_r + :math:`\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\theta` 1314 curl\_vp\_grad\_vp\_theta + :math:`\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\phi` 1315 curl\_vp\_grad\_vp\_phi + :math:`\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_r` 1316 curl\_vm\_grad\_vm\_r + :math:`\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\theta` 1317 curl\_vm\_grad\_vm\_theta + :math:`\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\phi` 1318 curl\_vm\_grad\_vm\_phi + :math:`\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\Theta\right)\right]_\theta` 1319 curl\_buoyancy\_force\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\Theta\right)\right]_\phi` 1320 curl\_buoyancy\_force\_phi + :math:`\left(\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\Theta\right)\right]_\theta\right)^2` 1321 curl\_buoyancy\_force\_theta\_squared + :math:`\left(\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\Theta\right)\right]_\phi\right)^2` 1322 curl\_buoyancy\_force\_phi\_squared + :math:`\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\overline{\Theta}\right)\right]_\theta` 1323 curl\_buoyancy\_pforce\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\overline{\Theta}\right)\right]_\phi` 1324 curl\_buoyancy\_pforce\_phi + :math:`\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\Theta'\right)\right]_\theta` 1325 curl\_buoyancy\_mforce\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\Theta'\right)\right]_\phi` 1326 curl\_buoyancy\_mforce\_phi + :math:`\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v})\right)\right]_r` 1327 curl\_coriolis\_force\_r + :math:`\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v})\right)\right]_\theta` 1328 curl\_coriolis\_force\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v})\right)\right]_\phi` 1329 curl\_coriolis\_force\_phi + :math:`\left(\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v})\right)\right]_r\right)^2` 1330 curl\_coriolis\_force\_r\_squared + :math:`\left(\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v})\right)\right]_\theta\right)^2` 1331 curl\_coriolis\_force\_theta\_squared + :math:`\left(\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v})\right)\right]_\phi\right)^2` 1332 curl\_coriolis\_force\_phi\_squared + :math:`\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v}')\right)\right]_r` 1333 curl\_coriolis\_pforce\_r + :math:`\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v}')\right)\right]_\theta` 1334 curl\_coriolis\_pforce\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\boldsymbol{v}')\right)\right]_\phi` 1335 curl\_coriolis\_pforce\_phi + :math:`\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\overline{\boldsymbol{\hat{v}}})\right)\right]_r` 1336 curl\_coriolis\_mforce\_r + :math:`\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\overline{\boldsymbol{\hat{v}}})\right)\right]_\theta` 1337 curl\_coriolis\_mforce\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(-c_1\mathrm{f}_1(\hat{\boldsymbol{z}}\times\overline{\boldsymbol{\hat{v}}})\right)\right]_\phi` 1338 curl\_coriolis\_mforce\_phi + :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_r` 1339 curl\_viscous\_force\_r + :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_\theta` 1340 curl\_viscous\_force\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_\phi` 1341 curl\_viscous\_force\_phi + :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_r` 1342 curl\_viscous\_pforce\_r + :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_\theta` 1343 curl\_viscous\_pforce\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_\phi` 1344 curl\_viscous\_pforce\_phi + :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_r` 1345 curl\_viscous\_mforce\_r + :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\theta` 1346 curl\_viscous\_mforce\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\phi` 1347 curl\_viscous\_mforce\_phi + :math:`\left\nabla\times\nabla P\right]_r` 1348 curl\_pressure\_force\_r + :math:`\left\nabla\times\nabla P\right]_\theta` 1349 curl\_pressure\_force\_theta + :math:`\left\nabla\times\nabla P\right]_\phi` 1350 curl\_pressure\_force\_phi + :math:`\left\nabla\times\nabla P\right]_r^2` 1351 curl\_pressure\_force\_r\_squared + :math:`\left\nabla\times\nabla P\right]_\theta^2` 1352 curl\_pressure\_force\_theta\_squared + :math:`\left\nabla\times\nabla P\right]_\phi^2` 1353 curl\_pressure\_force\_phi\_squared + :math:`\left\nabla\times\nabla P\right]_r` 1354 curl\_pressure\_pforce\_r + :math:`\left\nabla\times\nabla P\right]_\theta` 1355 curl\_pressure\_pforce\_theta + :math:`\left\nabla\times\nabla P\right]_\phi` 1356 curl\_pressure\_pforce\_phi + :math:`\left\nabla\times\nabla P\right]_r` 1357 curl\_pressure\_mforce\_r + :math:`\left\nabla\times\nabla P\right]_\theta` 1358 curl\_pressure\_mforce\_theta + :math:`\left\nabla\times\nabla P\right]_\phi` 1359 curl\_pressure\_mforce\_phi + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_r` 1360 curl\_j\_cross\_b\_r + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\theta` 1361 curl\_j\_cross\_b\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\phi` 1362 curl\_j\_cross\_b\_phi + :math:`\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_r\right)^2` 1363 curl\_j\_cross\_b\_r\_squared + :math:`\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\theta\right)^2` 1364 curl\_j\_cross\_b\_theta\_squared + :math:`\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\phi\right)^2` 1365 curl\_j\_cross\_b\_phi\_squared + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_r` 1366 curl\_jp\_cross\_bm\_r + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_\theta` 1367 curl\_jp\_cross\_bm\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_\phi` 1368 curl\_jp\_cross\_bm\_phi + :math:`\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_r` 1369 curl\_jm\_cross\_bp\_r + :math:`\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_\theta` 1370 curl\_jm\_cross\_bp\_theta + :math:`\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_\phi` 1371 curl\_jm\_cross\_bp\_phi + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_r` 1372 curl\_jm\_cross\_bm\_r + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_\theta` 1373 curl\_jm\_cross\_bm\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_\phi` 1374 curl\_jm\_cross\_bm\_phi + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_r` 1374 curl\_jp\_cross\_bp\_r + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_\theta` 1375 curl\_jp\_cross\_bp\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_\phi` 1376 curl\_jp\_cross\_bp\_phi +===================================================================================================================================================================== ====== ======================================= diff --git a/doc/source/diagnostic_codes/qcodes.rst b/doc/source/diagnostic_codes/qcodes.rst index cc3e71ed0..ad95c76af 100644 --- a/doc/source/diagnostic_codes/qcodes.rst +++ b/doc/source/diagnostic_codes/qcodes.rst @@ -17,6 +17,7 @@ Output Quantity Codes current_density magnetic_energy momentum_equation + curl_momentum_equation thermal_equation induction_equation amom_equation diff --git a/src/Diagnostics/curl_momentum_equation_codes.F b/src/Diagnostics/curl_momentum_equation_codes.F index 9cbf86669..17d4fc752 100644 --- a/src/Diagnostics/curl_momentum_equation_codes.F +++ b/src/Diagnostics/curl_momentum_equation_codes.F @@ -94,21 +94,21 @@ Integer, Parameter :: curl_viscous_mforce_phi = curl_mom_eq_off+47 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\phi$ ! Pressure forces - Integer, Parameter :: curl_pressure_force_r = curl_mom_eq_off+48 ! :tex: - Integer, Parameter :: curl_pressure_force_theta = curl_mom_eq_off+49 ! :tex: - Integer, Parameter :: curl_pressure_force_phi = curl_mom_eq_off+50 ! :tex: + Integer, Parameter :: curl_pressure_force_r = curl_mom_eq_off+48 ! :tex: $\left\nabla\times\nabla P\right]_r$ + Integer, Parameter :: curl_pressure_force_theta = curl_mom_eq_off+49 ! :tex: $\left\nabla\times\nabla P\right]_\theta$ + Integer, Parameter :: curl_pressure_force_phi = curl_mom_eq_off+50 ! :tex: $\left\nabla\times\nabla P\right]_\phi$ - Integer, Parameter :: curl_pressure_force_r_squared = curl_mom_eq_off+51 ! :tex: - Integer, Parameter :: curl_pressure_force_theta_squared = curl_mom_eq_off+52 ! :tex: - Integer, Parameter :: curl_pressure_force_phi_squared = curl_mom_eq_off+53 ! :tex: + Integer, Parameter :: curl_pressure_force_r_squared = curl_mom_eq_off+51 ! :tex: $\left\nabla\times\nabla P\right]_r^2$ + Integer, Parameter :: curl_pressure_force_theta_squared = curl_mom_eq_off+52 ! :tex: $\left\nabla\times\nabla P\right]_\theta^2$ + Integer, Parameter :: curl_pressure_force_phi_squared = curl_mom_eq_off+53 ! :tex: $\left\nabla\times\nabla P\right]_\phi^2$ -! Integer, Parameter :: curl_pressure_pforce_r = curl_mom_eq_off+54 ! :tex: -! Integer, Parameter :: curl_pressure_pforce_theta = curl_mom_eq_off+55 ! :tex: -! Integer, Parameter :: curl_pressure_pforce_phi = curl_mom_eq_off+56 ! :tex: +! Integer, Parameter :: curl_pressure_pforce_r = curl_mom_eq_off+54 ! :tex: $\left\nabla\times\nabla P\right]_r$ +! Integer, Parameter :: curl_pressure_pforce_theta = curl_mom_eq_off+55 ! :tex: $\left\nabla\times\nabla P\right]_\theta$ +! Integer, Parameter :: curl_pressure_pforce_phi = curl_mom_eq_off+56 ! :tex: $\left\nabla\times\nabla P\right]_\phi$ -! Integer, Parameter :: curl_pressure_mforce_r = curl_mom_eq_off+57 ! :tex: -! Integer, Parameter :: curl_pressure_mforce_theta = curl_mom_eq_off+58 ! :tex: -! Integer, Parameter :: curl_pressure_mforce_phi = curl_mom_eq_off+59 ! :tex: +! Integer, Parameter :: curl_pressure_mforce_r = curl_mom_eq_off+57 ! :tex: $\left\nabla\times\nabla P\right]_r$ +! Integer, Parameter :: curl_pressure_mforce_theta = curl_mom_eq_off+58 ! :tex: $\left\nabla\times\nabla P\right]_\theta$ +! Integer, Parameter :: curl_pressure_mforce_phi = curl_mom_eq_off+59 ! :tex: $\left\nabla\times\nabla P\right]_\phi$ !/////////////////////////// Lorentz forces /////////////////////////////// From 192d1ca6f8b389f369f7326a6ee1567e41b083ff Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Sun, 19 Jul 2026 17:28:21 +0000 Subject: [PATCH 08/13] Adding an untested example of turning the calculated gradients of the viscous force into a curl. --- src/Diagnostics/Diagnostics_Curl_Momentum.F90 | 37 +++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/src/Diagnostics/Diagnostics_Curl_Momentum.F90 b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 index e3d3d5681..5d77b3d5c 100644 --- a/src/Diagnostics/Diagnostics_Curl_Momentum.F90 +++ b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 @@ -38,6 +38,7 @@ Subroutine Compute_Curl_Momentum_Forces(buffer) Call Compute_Curl_Magnetic_Force(buffer) Call Compute_Curl_Coriolis_Force(buffer) Call Compute_Curl_Pressure_Force(buffer) + Call Compute_Curl_Viscous_Force(buffer) End Subroutine Subroutine Compute_Curl_Advection_Force(buffer) @@ -549,6 +550,42 @@ Subroutine Compute_Curl_Pressure_Force(buffer) End Subroutine Compute_Curl_Pressure_Force + Subroutine Compute_Curl_Viscous_Force(buffer) + Implicit None + Real*8, Intent(InOut) :: buffer(1:,my_r%min:,my_theta%min:,1:) + Integer :: r, k, t + + !!!!!!!!!!!!!!!!!!!!!!!!!!! Viscous Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + If (compute_quantity(curl_viscous_force_r)) Then + DO_PSI + qty(PSI) = One_Over_R(r)*(VFDBUFF(PSI,dvf_p_dt) + & + cottheta(t)*vforce_buffer(PSI,vfp_r) - & + csctheta(t)*VFDBUFF(PSI,dvf_t_dp)) + END_DO + Call Add_Quantity(qty) + Endif + + If (compute_quantity(curl_viscous_force_theta)) Then + DO_PSI + qty(PSI) = One_Over_R(r)*(csctheta(t)*VFDBUFF(PSI,dvf_r_dp) - & + vforce_buffer(PSI,vfp_p)) - & + VFDBUFF(PSI,dvf_p_dr) + END_DO + Call Add_Quantity(qty) + Endif + + If (compute_quantity(curl_viscous_force_phi)) Then + DO_PSI + qty(PSI) = VFDBUFF(PSI,dvf_t_dr) + & + One_Over_R(r)*(vforce_buffer(PSI,vfp_t) - & + VFDBUFF(PSI,dvf_r_dt)) + END_DO + Call Add_Quantity(qty) + Endif + + End Subroutine Compute_Curl_Viscous_Force + Subroutine Initialize_Grad_Viscous_Force() Implicit None integer :: nvfind, nvfdind, vfoff From ce218981e494ea67431b88f071ea21dfe82275b2 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Wed, 22 Jul 2026 13:43:24 -0400 Subject: [PATCH 09/13] Remove comments that have been resolved by #633 --- src/Diagnostics/Diagnostics_Curl_Momentum.F90 | 1 - src/Diagnostics/Diagnostics_Second_Derivatives.F90 | 1 - 2 files changed, 2 deletions(-) diff --git a/src/Diagnostics/Diagnostics_Curl_Momentum.F90 b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 index 5d77b3d5c..b63f7155a 100644 --- a/src/Diagnostics/Diagnostics_Curl_Momentum.F90 +++ b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 @@ -781,7 +781,6 @@ Subroutine Grad_Viscous_Force() endif enddo call gridcp%from_spectral(d_vforce_buffer%p1a, d_vforce_buffer%p1b) - ! FIXME: different to second derivatives d_vforce_buffer%p1a = d_vforce_buffer%p1b else do i = 1, nvffields diff --git a/src/Diagnostics/Diagnostics_Second_Derivatives.F90 b/src/Diagnostics/Diagnostics_Second_Derivatives.F90 index bc752f0cf..27ae6b8ec 100644 --- a/src/Diagnostics/Diagnostics_Second_Derivatives.F90 +++ b/src/Diagnostics/Diagnostics_Second_Derivatives.F90 @@ -400,7 +400,6 @@ Subroutine Compute_Second_Derivatives(inbuffer) Call d_by_dx3d3(i,j,d2buffer%p1a,1) Enddo Endif - ! CRW: Is this right in finite differences? d2buffer%p1a=d2buffer%p1b Call d2buffer%deconstruct('p1b') From 3d31fa07a30b7f8b027e953421ba2216e21eecb8 Mon Sep 17 00:00:00 2001 From: Ayesha Sarwar <2831664s@euclid-32.maths.gla.ac.uk> Date: Wed, 22 Jul 2026 19:42:56 +0100 Subject: [PATCH 10/13] Question about using new curl momentum diagnostics module Hi Cian, I have some questions: Firstly, I am unsure about the correct formulation of the curl of the pressure term. The current implementation does not appear to be mathematically correct The most important question is how I can use this new module now. Previously, I was simply adding the little routines directly in the Custom diagnostics, but with this new module, I am not sure what the correct procedure is for enabling and using it in the code. Thanks, Ayesha --- src/Diagnostics/Diagnostics_Curl_Momentum.F90 | 116 ++++++++++-------- .../curl_momentum_equation_codes.F | 60 ++++----- 2 files changed, 94 insertions(+), 82 deletions(-) diff --git a/src/Diagnostics/Diagnostics_Curl_Momentum.F90 b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 index b63f7155a..07a4103f2 100644 --- a/src/Diagnostics/Diagnostics_Curl_Momentum.F90 +++ b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 @@ -82,10 +82,9 @@ Subroutine Compute_Curl_Advection_Force(buffer) Call Add_Quantity(qty) Endif - ! New squared + sqrt quantity (i.e. |curl_v_grad_v_r|) If (compute_quantity(curl_v_grad_v_r_squared)) Then DO_PSI - qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + qty(PSI) = qty(PSI)*qty(PSI) END_DO Call Add_Quantity(qty) Endif @@ -127,7 +126,7 @@ Subroutine Compute_Curl_Advection_Force(buffer) Endif If (compute_quantity(curl_v_grad_v_theta_squared)) Then DO_PSI - qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + qty(PSI) = qty(PSI)*qty(PSI) END_DO Call Add_Quantity(qty) Endif @@ -172,7 +171,7 @@ Subroutine Compute_Curl_Advection_Force(buffer) Endif If (compute_quantity(curl_v_grad_v_phi_squared)) Then DO_PSI - qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + qty(PSI) = qty(PSI)*qty(PSI) END_DO Call Add_Quantity(qty) Endif @@ -198,10 +197,9 @@ Subroutine Compute_Curl_Buoyancy_Force(buffer) Call Add_Quantity(qty) Endif - ! New squared + sqrt quantity (i.e. |curl_buoyancy_force_theta|) If (compute_quantity(curl_buoyancy_force_theta_squared)) Then DO_PSI - qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + qty(PSI) = qty(PSI)*qty(PSI) END_DO Call Add_Quantity(qty) Endif @@ -216,10 +214,9 @@ Subroutine Compute_Curl_Buoyancy_Force(buffer) Call Add_Quantity(qty) Endif - ! New squared + sqrt quantity (i.e. |curl_buoyancy_force_phi|) If (compute_quantity(curl_buoyancy_force_phi_squared)) Then DO_PSI - qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + qty(PSI) = qty(PSI)*qty(PSI) END_DO Call Add_Quantity(qty) Endif @@ -287,10 +284,9 @@ Subroutine Compute_Curl_Magnetic_Force(buffer) Call Add_Quantity(qty) Endif - ! New squared + sqrt quantity (i.e. |curl_j_cross_b_r|) If (compute_quantity(curl_j_cross_b_r_squared)) Then DO_PSI - qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + qty(PSI) = qty(PSI)*qty(PSI) END_DO Call Add_Quantity(qty) Endif @@ -343,10 +339,9 @@ Subroutine Compute_Curl_Magnetic_Force(buffer) Call Add_Quantity(qty) Endif - ! New squared + sqrt quantity (i.e. |curl_j_cross_b_theta|) If (compute_quantity(curl_j_cross_b_theta_squared)) Then DO_PSI - qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + qty(PSI) = qty(PSI)*qty(PSI) END_DO Call Add_Quantity(qty) Endif @@ -396,10 +391,9 @@ Subroutine Compute_Curl_Magnetic_Force(buffer) Call Add_Quantity(qty) Endif - ! New squared + sqrt quantity (i.e. |curl_j_cross_b_phi|) If (compute_quantity(curl_j_cross_b_phi_squared)) Then DO_PSI - qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + qty(PSI) = qty(PSI)*qty(PSI) END_DO Call Add_Quantity(qty) Endif @@ -426,10 +420,9 @@ Subroutine Compute_Curl_Coriolis_Force(buffer) Call Add_Quantity(qty) Endif - ! New squared + sqrt quantity (i.e. |curl_coriolis_force_r|) If (compute_quantity(curl_coriolis_force_r_squared)) Then DO_PSI - qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + qty(PSI) = qty(PSI)*qty(PSI) END_DO Call Add_Quantity(qty) Endif @@ -447,10 +440,9 @@ Subroutine Compute_Curl_Coriolis_Force(buffer) Call Add_Quantity(qty) Endif - ! New squared + sqrt quantity (i.e. |curl_coriolis_force_theta|) If (compute_quantity(curl_coriolis_force_theta_squared)) Then DO_PSI - qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + qty(PSI) = qty(PSI)*qty(PSI) END_DO Call Add_Quantity(qty) Endif @@ -467,10 +459,9 @@ Subroutine Compute_Curl_Coriolis_Force(buffer) Call Add_Quantity(qty) Endif - ! New squared + sqrt quantity (i.e. |curl_coriolis_force_phi|) If (compute_quantity(curl_coriolis_force_phi_squared)) Then DO_PSI - qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + qty(PSI) = qty(PSI)*qty(PSI) END_DO Call Add_Quantity(qty) Endif @@ -499,10 +490,9 @@ Subroutine Compute_Curl_Pressure_Force(buffer) Call Add_Quantity(qty) Endif - ! New squared + sqrt quantity (i.e. |curl_viscous_force_r|) If (compute_quantity(curl_pressure_force_r_squared)) Then DO_PSI - qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + qty(PSI) = qty(PSI)*qty(PSI) END_DO Call Add_Quantity(qty) Endif @@ -518,10 +508,9 @@ Subroutine Compute_Curl_Pressure_Force(buffer) Call Add_Quantity(qty) Endif - ! New squared + sqrt quantity (i.e. |curl_pressure_force_theta|) If (compute_quantity(curl_pressure_force_theta_squared)) Then DO_PSI - qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + qty(PSI) = qty(PSI)*qty(PSI) END_DO Call Add_Quantity(qty) Endif @@ -539,10 +528,9 @@ Subroutine Compute_Curl_Pressure_Force(buffer) Call Add_Quantity(qty) Endif - ! New squared + sqrt quantity (i.e. |curl_viscous_force_phi|) If (compute_quantity(curl_pressure_force_phi_squared)) Then DO_PSI - qty(PSI) = sqrt(qty(PSI)*qty(PSI)) + qty(PSI) = qty(PSI)*qty(PSI) END_DO Call Add_Quantity(qty) Endif @@ -556,33 +544,57 @@ Subroutine Compute_Curl_Viscous_Force(buffer) Integer :: r, k, t !!!!!!!!!!!!!!!!!!!!!!!!!!! Viscous Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + If (compute_quantity(curl_viscous_force_r) .or. compute_quantity(curl_viscous_force_r_squared)) Then + If (compute_quantity(curl_viscous_force_r)) Then + DO_PSI + qty(PSI) = One_Over_R(r)*(VFDBUFF(PSI,dvf_p_dt) + & + cottheta(t)*vforce_buffer(PSI,vfp_r) - & + csctheta(t)*VFDBUFF(PSI,dvf_t_dp)) + END_DO + Call Add_Quantity(qty) + Endif + If (compute_quantity(curl_viscous_force_r_squared)) Then + DO_PSI + qty(PSI) = qty(PSI)*qty(PSI) + END_DO + Call Add_Quantity(qty) + Endif + Endif + + If (compute_quantity(curl_viscous_force_theta) .or. compute_quantity(curl_viscous_force_theta_squared)) Then + If (compute_quantity(curl_viscous_force_theta)) Then + DO_PSI + qty(PSI) = One_Over_R(r)*(csctheta(t)*VFDBUFF(PSI,dvf_r_dp) - & + vforce_buffer(PSI,vfp_p)) - & + VFDBUFF(PSI,dvf_p_dr) + END_DO + Call Add_Quantity(qty) + Endif + If (compute_quantity(curl_viscous_force_theta_squared)) Then + DO_PSI + qty(PSI) = qty(PSI)*qty(PSI) + END_DO + Call Add_Quantity(qty) + Endif + Endif - If (compute_quantity(curl_viscous_force_r)) Then - DO_PSI - qty(PSI) = One_Over_R(r)*(VFDBUFF(PSI,dvf_p_dt) + & - cottheta(t)*vforce_buffer(PSI,vfp_r) - & - csctheta(t)*VFDBUFF(PSI,dvf_t_dp)) - END_DO - Call Add_Quantity(qty) - Endif - - If (compute_quantity(curl_viscous_force_theta)) Then - DO_PSI - qty(PSI) = One_Over_R(r)*(csctheta(t)*VFDBUFF(PSI,dvf_r_dp) - & - vforce_buffer(PSI,vfp_p)) - & - VFDBUFF(PSI,dvf_p_dr) - END_DO - Call Add_Quantity(qty) - Endif - - If (compute_quantity(curl_viscous_force_phi)) Then - DO_PSI - qty(PSI) = VFDBUFF(PSI,dvf_t_dr) + & - One_Over_R(r)*(vforce_buffer(PSI,vfp_t) - & - VFDBUFF(PSI,dvf_r_dt)) - END_DO - Call Add_Quantity(qty) - Endif + If (compute_quantity(curl_viscous_force_phi) .or. compute_quantity(curl_viscous_force_phi_squared)) Then + If (compute_quantity(curl_viscous_force_phi)) Then + DO_PSI + qty(PSI) = VFDBUFF(PSI,dvf_t_dr) + & + One_Over_R(r)*(vforce_buffer(PSI,vfp_t) - & + VFDBUFF(PSI,dvf_r_dt)) + END_DO + Call Add_Quantity(qty) + Endif + If (compute_quantity(curl_viscous_force_phi_squared)) Then + DO_PSI + qty(PSI) = qty(PSI)*qty(PSI) + END_DO + Call Add_Quantity(qty) + Endif + Endif End Subroutine Compute_Curl_Viscous_Force diff --git a/src/Diagnostics/curl_momentum_equation_codes.F b/src/Diagnostics/curl_momentum_equation_codes.F index 17d4fc752..d75f12bf0 100644 --- a/src/Diagnostics/curl_momentum_equation_codes.F +++ b/src/Diagnostics/curl_momentum_equation_codes.F @@ -23,29 +23,29 @@ ! NOTE: ADVECTION TERMS ARE SCALED BY DENSITY (so that they represent a force density) Integer, Parameter :: curl_mom_eq_off = mom_eq_off+100 ! :OFFSET CODE: - Integer, Parameter :: curl_v_grad_v_r = curl_mom_eq_off+1 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_r$ - Integer, Parameter :: curl_v_grad_v_theta = curl_mom_eq_off+2 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\theta$ - Integer, Parameter :: curl_v_grad_v_phi = curl_mom_eq_off+3 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\phi$ + Integer, Parameter :: curl_v_grad_v_r = curl_mom_eq_off+1 ! :tex: $\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]\right]_r$ + Integer, Parameter :: curl_v_grad_v_theta = curl_mom_eq_off+2 ! :tex: $\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]\right]_\theta$ + Integer, Parameter :: curl_v_grad_v_phi = curl_mom_eq_off+3 ! :tex: $\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]\right]_\phi$ - Integer, Parameter :: curl_v_grad_v_r_squared = curl_mom_eq_off+4 ! :tex: $\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_r\right)^2$ - Integer, Parameter :: curl_v_grad_v_theta_squared = curl_mom_eq_off+5 ! :tex: $\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\theta\right)^2$ - Integer, Parameter :: curl_v_grad_v_phi_squared = curl_mom_eq_off+6 ! :tex: $\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\phi\right)^2$ + Integer, Parameter :: curl_v_grad_v_r_squared = curl_mom_eq_off+4 ! :tex: $\left(\left[\nabla\times\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]\right)\right]_r\right)^2$ + Integer, Parameter :: curl_v_grad_v_theta_squared = curl_mom_eq_off+5 ! :tex: $\left(\left[\nabla\times\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]\right)\right]_\theta\right)^2$ + Integer, Parameter :: curl_v_grad_v_phi_squared = curl_mom_eq_off+6 ! :tex: $\left(\left[\nabla\times\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]\right)\right]_\phi\right)^2$ -! Integer, Parameter :: curl_vp_grad_vm_r = curl_mom_eq_off+7 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_r$ -! Integer, Parameter :: curl_vp_grad_vm_theta = curl_mom_eq_off+8 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\theta$ -! Integer, Parameter :: curl_vp_grad_vm_phi = curl_mom_eq_off+9 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\phi$ +! Integer, Parameter :: curl_vp_grad_vm_r = curl_mom_eq_off+7 ! :tex: $\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]\right]_r$ +! Integer, Parameter :: curl_vp_grad_vm_theta = curl_mom_eq_off+8 ! :tex: $\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]\right]_\theta$ +! Integer, Parameter :: curl_vp_grad_vm_phi = curl_mom_eq_off+9 ! :tex: $\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]\right]_\phi$ -! Integer, Parameter :: curl_vm_grad_vp_r = curl_mom_eq_off+10 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_r$ -! Integer, Parameter :: curl_vm_grad_vp_theta = curl_mom_eq_off+11 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\theta$ -! Integer, Parameter :: curl_vm_grad_vp_phi = curl_mom_eq_off+12 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\phi$ +! Integer, Parameter :: curl_vm_grad_vp_r = curl_mom_eq_off+10 ! :tex: $\left[\nabla\times\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]\right]_r$ +! Integer, Parameter :: curl_vm_grad_vp_theta = curl_mom_eq_off+11 ! :tex: $\left[\nabla\times\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]\right]_\theta$ +! Integer, Parameter :: curl_vm_grad_vp_phi = curl_mom_eq_off+12 ! :tex: $\left[\nabla\times\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]\right]_\phi$ -! Integer, Parameter :: curl_vp_grad_vp_r = curl_mom_eq_off+13 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_r$ -! Integer, Parameter :: curl_vp_grad_vp_theta = curl_mom_eq_off+14 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\theta$ -! Integer, Parameter :: curl_vp_grad_vp_phi = curl_mom_eq_off+15 ! :tex: $\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\phi$ +! Integer, Parameter :: curl_vp_grad_vp_r = curl_mom_eq_off+13 ! :tex: $\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]\right]_r$ +! Integer, Parameter :: curl_vp_grad_vp_theta = curl_mom_eq_off+14 ! :tex: $\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]\right]_\theta$ +! Integer, Parameter :: curl_vp_grad_vp_phi = curl_mom_eq_off+15 ! :tex: $\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]\right]_\phi$ -! Integer, Parameter :: curl_vm_grad_vm_r = curl_mom_eq_off+16 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_r$ -! Integer, Parameter :: curl_vm_grad_vm_theta = curl_mom_eq_off+17 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\theta$ -! Integer, Parameter :: curl_vm_grad_vm_phi = curl_mom_eq_off+18 ! :tex: $\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\phi$ +! Integer, Parameter :: curl_vm_grad_vm_r = curl_mom_eq_off+16 ! :tex: $\left[\nabla\times\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]\right]_r$ +! Integer, Parameter :: curl_vm_grad_vm_theta = curl_mom_eq_off+17 ! :tex: $\left[\nabla\times\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]\right]_\theta$ +! Integer, Parameter :: curl_vm_grad_vm_phi = curl_mom_eq_off+18 ! :tex: $\left[\nabla\times\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]\right]_\phi$ !///////////////////////////////////////////////////////////// ! Linear forces @@ -94,21 +94,21 @@ Integer, Parameter :: curl_viscous_mforce_phi = curl_mom_eq_off+47 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\phi$ ! Pressure forces - Integer, Parameter :: curl_pressure_force_r = curl_mom_eq_off+48 ! :tex: $\left\nabla\times\nabla P\right]_r$ - Integer, Parameter :: curl_pressure_force_theta = curl_mom_eq_off+49 ! :tex: $\left\nabla\times\nabla P\right]_\theta$ - Integer, Parameter :: curl_pressure_force_phi = curl_mom_eq_off+50 ! :tex: $\left\nabla\times\nabla P\right]_\phi$ + Integer, Parameter :: curl_pressure_force_r = curl_mom_eq_off+48 ! :tex: $\left[\nabla\times\nabla P\right]_r$ + Integer, Parameter :: curl_pressure_force_theta = curl_mom_eq_off+49 ! :tex: $\left[\nabla\times\nabla P\right]_\theta$ + Integer, Parameter :: curl_pressure_force_phi = curl_mom_eq_off+50 ! :tex: $\left[\nabla\times\nabla P\right]_\phi$ - Integer, Parameter :: curl_pressure_force_r_squared = curl_mom_eq_off+51 ! :tex: $\left\nabla\times\nabla P\right]_r^2$ - Integer, Parameter :: curl_pressure_force_theta_squared = curl_mom_eq_off+52 ! :tex: $\left\nabla\times\nabla P\right]_\theta^2$ - Integer, Parameter :: curl_pressure_force_phi_squared = curl_mom_eq_off+53 ! :tex: $\left\nabla\times\nabla P\right]_\phi^2$ + Integer, Parameter :: curl_pressure_force_r_squared = curl_mom_eq_off+51 ! :tex: $\left[\nabla\times\nabla P\right]_r^2$ + Integer, Parameter :: curl_pressure_force_theta_squared = curl_mom_eq_off+52 ! :tex: $\left[\nabla\times\nabla P\right]_\theta^2$ + Integer, Parameter :: curl_pressure_force_phi_squared = curl_mom_eq_off+53 ! :tex: $\left[\nabla\times\nabla P\right]_\phi^2$ -! Integer, Parameter :: curl_pressure_pforce_r = curl_mom_eq_off+54 ! :tex: $\left\nabla\times\nabla P\right]_r$ -! Integer, Parameter :: curl_pressure_pforce_theta = curl_mom_eq_off+55 ! :tex: $\left\nabla\times\nabla P\right]_\theta$ -! Integer, Parameter :: curl_pressure_pforce_phi = curl_mom_eq_off+56 ! :tex: $\left\nabla\times\nabla P\right]_\phi$ +! Integer, Parameter :: curl_pressure_pforce_r = curl_mom_eq_off+54 ! :tex: $\left[\nabla\times\nabla P\right]_r$ +! Integer, Parameter :: curl_pressure_pforce_theta = curl_mom_eq_off+55 ! :tex: $\left[\nabla\times\nabla P\right]_\theta$ +! Integer, Parameter :: curl_pressure_pforce_phi = curl_mom_eq_off+56 ! :tex: $\left[\nabla\times\nabla P\right]_\phi$ -! Integer, Parameter :: curl_pressure_mforce_r = curl_mom_eq_off+57 ! :tex: $\left\nabla\times\nabla P\right]_r$ -! Integer, Parameter :: curl_pressure_mforce_theta = curl_mom_eq_off+58 ! :tex: $\left\nabla\times\nabla P\right]_\theta$ -! Integer, Parameter :: curl_pressure_mforce_phi = curl_mom_eq_off+59 ! :tex: $\left\nabla\times\nabla P\right]_\phi$ +! Integer, Parameter :: curl_pressure_mforce_r = curl_mom_eq_off+57 ! :tex: $\left[\nabla\times\nabla P\right]_r$ +! Integer, Parameter :: curl_pressure_mforce_theta = curl_mom_eq_off+58 ! :tex: $\left[\nabla\times\nabla P\right]_\theta$ +! Integer, Parameter :: curl_pressure_mforce_phi = curl_mom_eq_off+59 ! :tex: $\left[\nabla\times\nabla P\right]_\phi$ !/////////////////////////// Lorentz forces /////////////////////////////// From 5382d47505cb274da7c4ced5eb149fd1edaa706a Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Wed, 22 Jul 2026 18:29:33 -0400 Subject: [PATCH 11/13] Fix for the failing build test. Needed to add the quantity codes to the header file. --- .../curl_momentum_equation.rst | 111 +++++++++--------- .../curl_momentum_equation_codes.F | 76 ++++++------ 2 files changed, 97 insertions(+), 90 deletions(-) diff --git a/doc/source/diagnostic_codes/curl_momentum_equation.rst b/doc/source/diagnostic_codes/curl_momentum_equation.rst index 0729d2762..64c4c53e3 100644 --- a/doc/source/diagnostic_codes/curl_momentum_equation.rst +++ b/doc/source/diagnostic_codes/curl_momentum_equation.rst @@ -2,24 +2,24 @@ Curl Momentum Equation ==================================================================== ===================================================================================================================================================================== ====== ======================================= - :math:`\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_r` 1301 curl\_v\_grad\_v\_r - :math:`\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\theta` 1302 curl\_v\_grad\_v\_theta - :math:`\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\phi` 1303 curl\_v\_grad\_v\_phi - :math:`\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_r\right)^2` 1304 curl\_v\_grad\_v\_r\_squared - :math:`\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\theta\right)^2` 1305 curl\_v\_grad\_v\_theta\_squared - :math:`\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]_\phi\right)^2` 1306 curl\_v\_grad\_v\_phi\_squared - :math:`\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_r` 1307 curl\_vp\_grad\_vm\_r - :math:`\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\theta` 1308 curl\_vp\_grad\_vm\_theta - :math:`\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\phi` 1309 curl\_vp\_grad\_vm\_phi - :math:`\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_r` 1310 curl\_vm\_grad\_vp\_r - :math:`\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\theta` 1311 curl\_vm\_grad\_vp\_theta - :math:`\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\phi` 1312 curl\_vm\_grad\_vp\_phi - :math:`\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_r` 1313 curl\_vp\_grad\_vp\_r - :math:`\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\theta` 1314 curl\_vp\_grad\_vp\_theta - :math:`\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]_\phi` 1315 curl\_vp\_grad\_vp\_phi - :math:`\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_r` 1316 curl\_vm\_grad\_vm\_r - :math:`\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\theta` 1317 curl\_vm\_grad\_vm\_theta - :math:`\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]_\phi` 1318 curl\_vm\_grad\_vm\_phi + :math:`\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]\right]_r` 1301 curl\_v\_grad\_v\_r + :math:`\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]\right]_\theta` 1302 curl\_v\_grad\_v\_theta + :math:`\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]\right]_\phi` 1303 curl\_v\_grad\_v\_phi + :math:`\left(\left[\nabla\times\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]\right)\right]_r\right)^2` 1304 curl\_v\_grad\_v\_r\_squared + :math:`\left(\left[\nabla\times\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]\right)\right]_\theta\right)^2` 1305 curl\_v\_grad\_v\_theta\_squared + :math:`\left(\left[\nabla\times\left(\mathrm{f}_1\left[\boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}\right]\right)\right]_\phi\right)^2` 1306 curl\_v\_grad\_v\_phi\_squared + :math:`\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]\right]_r` 1307 curl\_vp\_grad\_vm\_r + :math:`\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]\right]_\theta` 1308 curl\_vp\_grad\_vm\_theta + :math:`\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]\right]_\phi` 1309 curl\_vp\_grad\_vm\_phi + :math:`\left[\nabla\times\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]\right]_r` 1310 curl\_vm\_grad\_vp\_r + :math:`\left[\nabla\times\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]\right]_\theta` 1311 curl\_vm\_grad\_vp\_theta + :math:`\left[\nabla\times\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]\right]_\phi` 1312 curl\_vm\_grad\_vp\_phi + :math:`\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]\right]_r` 1313 curl\_vp\_grad\_vp\_r + :math:`\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]\right]_\theta` 1314 curl\_vp\_grad\_vp\_theta + :math:`\left[\nabla\times\mathrm{f}_1\left[\boldsymbol{v'}\cdot\boldsymbol{\nabla}\boldsymbol{v'}\right]\right]_\phi` 1315 curl\_vp\_grad\_vp\_phi + :math:`\left[\nabla\times\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]\right]_r` 1316 curl\_vm\_grad\_vm\_r + :math:`\left[\nabla\times\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]\right]_\theta` 1317 curl\_vm\_grad\_vm\_theta + :math:`\left[\nabla\times\mathrm{f}_1\left[\overline{\boldsymbol{v}}\cdot\boldsymbol{\nabla}\overline{\boldsymbol{v}}\right]\right]_\phi` 1318 curl\_vm\_grad\_vm\_phi :math:`\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\Theta\right)\right]_\theta` 1319 curl\_buoyancy\_force\_theta :math:`\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\Theta\right)\right]_\phi` 1320 curl\_buoyancy\_force\_phi :math:`\left(\left[\boldsymbol{\nabla}\times\left(c_2\mathrm{f}_2\Theta\right)\right]_\theta\right)^2` 1321 curl\_buoyancy\_force\_theta\_squared @@ -43,40 +43,43 @@ Curl Momentum Equation :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_r` 1339 curl\_viscous\_force\_r :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_\theta` 1340 curl\_viscous\_force\_theta :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_\phi` 1341 curl\_viscous\_force\_phi - :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_r` 1342 curl\_viscous\_pforce\_r - :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_\theta` 1343 curl\_viscous\_pforce\_theta - :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_\phi` 1344 curl\_viscous\_pforce\_phi - :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_r` 1345 curl\_viscous\_mforce\_r - :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\theta` 1346 curl\_viscous\_mforce\_theta - :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\phi` 1347 curl\_viscous\_mforce\_phi - :math:`\left\nabla\times\nabla P\right]_r` 1348 curl\_pressure\_force\_r - :math:`\left\nabla\times\nabla P\right]_\theta` 1349 curl\_pressure\_force\_theta - :math:`\left\nabla\times\nabla P\right]_\phi` 1350 curl\_pressure\_force\_phi - :math:`\left\nabla\times\nabla P\right]_r^2` 1351 curl\_pressure\_force\_r\_squared - :math:`\left\nabla\times\nabla P\right]_\theta^2` 1352 curl\_pressure\_force\_theta\_squared - :math:`\left\nabla\times\nabla P\right]_\phi^2` 1353 curl\_pressure\_force\_phi\_squared - :math:`\left\nabla\times\nabla P\right]_r` 1354 curl\_pressure\_pforce\_r - :math:`\left\nabla\times\nabla P\right]_\theta` 1355 curl\_pressure\_pforce\_theta - :math:`\left\nabla\times\nabla P\right]_\phi` 1356 curl\_pressure\_pforce\_phi - :math:`\left\nabla\times\nabla P\right]_r` 1357 curl\_pressure\_mforce\_r - :math:`\left\nabla\times\nabla P\right]_\theta` 1358 curl\_pressure\_mforce\_theta - :math:`\left\nabla\times\nabla P\right]_\phi` 1359 curl\_pressure\_mforce\_phi - :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_r` 1360 curl\_j\_cross\_b\_r - :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\theta` 1361 curl\_j\_cross\_b\_theta - :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\phi` 1362 curl\_j\_cross\_b\_phi - :math:`\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_r\right)^2` 1363 curl\_j\_cross\_b\_r\_squared - :math:`\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\theta\right)^2` 1364 curl\_j\_cross\_b\_theta\_squared - :math:`\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\phi\right)^2` 1365 curl\_j\_cross\_b\_phi\_squared - :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_r` 1366 curl\_jp\_cross\_bm\_r - :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_\theta` 1367 curl\_jp\_cross\_bm\_theta - :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_\phi` 1368 curl\_jp\_cross\_bm\_phi - :math:`\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_r` 1369 curl\_jm\_cross\_bp\_r - :math:`\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_\theta` 1370 curl\_jm\_cross\_bp\_theta - :math:`\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_\phi` 1371 curl\_jm\_cross\_bp\_phi - :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_r` 1372 curl\_jm\_cross\_bm\_r - :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_\theta` 1373 curl\_jm\_cross\_bm\_theta - :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_\phi` 1374 curl\_jm\_cross\_bm\_phi - :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_r` 1374 curl\_jp\_cross\_bp\_r - :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_\theta` 1375 curl\_jp\_cross\_bp\_theta - :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_\phi` 1376 curl\_jp\_cross\_bp\_phi + :math:`\left(\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_r\right)^2` 1342 curl\_viscous\_force\_r\_squared + :math:`\left(\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_\theta\right)^2` 1343 curl\_viscous\_force\_theta\_squared + :math:`\left(\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_\phi\right)^2` 1344 curl\_viscous\_force\_phi\_squared + :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_r` 1345 curl\_viscous\_pforce\_r + :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_\theta` 1346 curl\_viscous\_pforce\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_\phi` 1347 curl\_viscous\_pforce\_phi + :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_r` 1351 curl\_viscous\_mforce\_r + :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\theta` 1352 curl\_viscous\_mforce\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\phi` 1353 curl\_viscous\_mforce\_phi + :math:`\left[\nabla\times\nabla P\right]_r` 1357 curl\_pressure\_force\_r + :math:`\left[\nabla\times\nabla P\right]_\theta` 1358 curl\_pressure\_force\_theta + :math:`\left[\nabla\times\nabla P\right]_\phi` 1359 curl\_pressure\_force\_phi + :math:`\left[\nabla\times\nabla P\right]_r^2` 1360 curl\_pressure\_force\_r\_squared + :math:`\left[\nabla\times\nabla P\right]_\theta^2` 1361 curl\_pressure\_force\_theta\_squared + :math:`\left[\nabla\times\nabla P\right]_\phi^2` 1362 curl\_pressure\_force\_phi\_squared + :math:`\left[\nabla\times\nabla P\right]_r` 1363 curl\_pressure\_pforce\_r + :math:`\left[\nabla\times\nabla P\right]_\theta` 1364 curl\_pressure\_pforce\_theta + :math:`\left[\nabla\times\nabla P\right]_\phi` 1365 curl\_pressure\_pforce\_phi + :math:`\left[\nabla\times\nabla P\right]_r` 1366 curl\_pressure\_mforce\_r + :math:`\left[\nabla\times\nabla P\right]_\theta` 1367 curl\_pressure\_mforce\_theta + :math:`\left[\nabla\times\nabla P\right]_\phi` 1368 curl\_pressure\_mforce\_phi + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_r` 1369 curl\_j\_cross\_b\_r + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\theta` 1370 curl\_j\_cross\_b\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\phi` 1371 curl\_j\_cross\_b\_phi + :math:`\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_r\right)^2` 1372 curl\_j\_cross\_b\_r\_squared + :math:`\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\theta\right)^2` 1373 curl\_j\_cross\_b\_theta\_squared + :math:`\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\phi\right)^2` 1374 curl\_j\_cross\_b\_phi\_squared + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_r` 1375 curl\_jp\_cross\_bm\_r + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_\theta` 1376 curl\_jp\_cross\_bm\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_\phi` 1377 curl\_jp\_cross\_bm\_phi + :math:`\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_r` 1378 curl\_jm\_cross\_bp\_r + :math:`\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_\theta` 1379 curl\_jm\_cross\_bp\_theta + :math:`\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_\phi` 1380 curl\_jm\_cross\_bp\_phi + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_r` 1381 curl\_jm\_cross\_bm\_r + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_\theta` 1382 curl\_jm\_cross\_bm\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_\phi` 1383 curl\_jm\_cross\_bm\_phi + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_r` 1384 curl\_jp\_cross\_bp\_r + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_\theta` 1385 curl\_jp\_cross\_bp\_theta + :math:`\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_\phi` 1386 curl\_jp\_cross\_bp\_phi ===================================================================================================================================================================== ====== ======================================= diff --git a/src/Diagnostics/curl_momentum_equation_codes.F b/src/Diagnostics/curl_momentum_equation_codes.F index d75f12bf0..d8452215f 100644 --- a/src/Diagnostics/curl_momentum_equation_codes.F +++ b/src/Diagnostics/curl_momentum_equation_codes.F @@ -85,30 +85,34 @@ Integer, Parameter :: curl_viscous_force_theta = curl_mom_eq_off+40 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_\theta$ Integer, Parameter :: curl_viscous_force_phi = curl_mom_eq_off+41 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_\phi$ - Integer, Parameter :: curl_viscous_pforce_r = curl_mom_eq_off+42 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_r$ - Integer, Parameter :: curl_viscous_pforce_theta = curl_mom_eq_off+43 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_\theta$ - Integer, Parameter :: curl_viscous_pforce_phi = curl_mom_eq_off+44 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_\phi$ + Integer, Parameter :: curl_viscous_force_r_squared = curl_mom_eq_off+42 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_r\right)^2$ + Integer, Parameter :: curl_viscous_force_theta_squared = curl_mom_eq_off+43 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_\theta\right)^2$ + Integer, Parameter :: curl_viscous_force_phi_squared = curl_mom_eq_off+44 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D})\right)\right]_\phi\right)^2$ - Integer, Parameter :: curl_viscous_mforce_r = curl_mom_eq_off+45 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_r$ - Integer, Parameter :: curl_viscous_mforce_theta = curl_mom_eq_off+46 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\theta$ - Integer, Parameter :: curl_viscous_mforce_phi = curl_mom_eq_off+47 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\phi$ + Integer, Parameter :: curl_viscous_pforce_r = curl_mom_eq_off+45 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_r$ + Integer, Parameter :: curl_viscous_pforce_theta = curl_mom_eq_off+46 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_\theta$ + Integer, Parameter :: curl_viscous_pforce_phi = curl_mom_eq_off+47 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal D'})\right)\right]_\phi$ + + Integer, Parameter :: curl_viscous_mforce_r = curl_mom_eq_off+51 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_r$ + Integer, Parameter :: curl_viscous_mforce_theta = curl_mom_eq_off+52 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\theta$ + Integer, Parameter :: curl_viscous_mforce_phi = curl_mom_eq_off+53 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_5\,(\boldsymbol{\nabla}\cdot\overline{\boldsymbol{\mathcal D}})\right)\right]_\phi$ ! Pressure forces - Integer, Parameter :: curl_pressure_force_r = curl_mom_eq_off+48 ! :tex: $\left[\nabla\times\nabla P\right]_r$ - Integer, Parameter :: curl_pressure_force_theta = curl_mom_eq_off+49 ! :tex: $\left[\nabla\times\nabla P\right]_\theta$ - Integer, Parameter :: curl_pressure_force_phi = curl_mom_eq_off+50 ! :tex: $\left[\nabla\times\nabla P\right]_\phi$ + Integer, Parameter :: curl_pressure_force_r = curl_mom_eq_off+57 ! :tex: $\left[\nabla\times\nabla P\right]_r$ + Integer, Parameter :: curl_pressure_force_theta = curl_mom_eq_off+58 ! :tex: $\left[\nabla\times\nabla P\right]_\theta$ + Integer, Parameter :: curl_pressure_force_phi = curl_mom_eq_off+59 ! :tex: $\left[\nabla\times\nabla P\right]_\phi$ - Integer, Parameter :: curl_pressure_force_r_squared = curl_mom_eq_off+51 ! :tex: $\left[\nabla\times\nabla P\right]_r^2$ - Integer, Parameter :: curl_pressure_force_theta_squared = curl_mom_eq_off+52 ! :tex: $\left[\nabla\times\nabla P\right]_\theta^2$ - Integer, Parameter :: curl_pressure_force_phi_squared = curl_mom_eq_off+53 ! :tex: $\left[\nabla\times\nabla P\right]_\phi^2$ + Integer, Parameter :: curl_pressure_force_r_squared = curl_mom_eq_off+60 ! :tex: $\left[\nabla\times\nabla P\right]_r^2$ + Integer, Parameter :: curl_pressure_force_theta_squared = curl_mom_eq_off+61 ! :tex: $\left[\nabla\times\nabla P\right]_\theta^2$ + Integer, Parameter :: curl_pressure_force_phi_squared = curl_mom_eq_off+62 ! :tex: $\left[\nabla\times\nabla P\right]_\phi^2$ -! Integer, Parameter :: curl_pressure_pforce_r = curl_mom_eq_off+54 ! :tex: $\left[\nabla\times\nabla P\right]_r$ -! Integer, Parameter :: curl_pressure_pforce_theta = curl_mom_eq_off+55 ! :tex: $\left[\nabla\times\nabla P\right]_\theta$ -! Integer, Parameter :: curl_pressure_pforce_phi = curl_mom_eq_off+56 ! :tex: $\left[\nabla\times\nabla P\right]_\phi$ +! Integer, Parameter :: curl_pressure_pforce_r = curl_mom_eq_off+63 ! :tex: $\left[\nabla\times\nabla P\right]_r$ +! Integer, Parameter :: curl_pressure_pforce_theta = curl_mom_eq_off+64 ! :tex: $\left[\nabla\times\nabla P\right]_\theta$ +! Integer, Parameter :: curl_pressure_pforce_phi = curl_mom_eq_off+65 ! :tex: $\left[\nabla\times\nabla P\right]_\phi$ -! Integer, Parameter :: curl_pressure_mforce_r = curl_mom_eq_off+57 ! :tex: $\left[\nabla\times\nabla P\right]_r$ -! Integer, Parameter :: curl_pressure_mforce_theta = curl_mom_eq_off+58 ! :tex: $\left[\nabla\times\nabla P\right]_\theta$ -! Integer, Parameter :: curl_pressure_mforce_phi = curl_mom_eq_off+59 ! :tex: $\left[\nabla\times\nabla P\right]_\phi$ +! Integer, Parameter :: curl_pressure_mforce_r = curl_mom_eq_off+66 ! :tex: $\left[\nabla\times\nabla P\right]_r$ +! Integer, Parameter :: curl_pressure_mforce_theta = curl_mom_eq_off+67 ! :tex: $\left[\nabla\times\nabla P\right]_\theta$ +! Integer, Parameter :: curl_pressure_mforce_phi = curl_mom_eq_off+68 ! :tex: $\left[\nabla\times\nabla P\right]_\phi$ !/////////////////////////// Lorentz forces /////////////////////////////// @@ -116,27 +120,27 @@ ! ref%Lorentz_Coeff = 1/4pi when dimensional, Pr/(Pr_m E) when nondimesional ! j (below) is shorthand for ref%Lorentz_Coeff*delxB (not quite the current density) - Integer, Parameter :: curl_j_cross_b_r = curl_mom_eq_off+60 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_r$ - Integer, Parameter :: curl_j_cross_b_theta = curl_mom_eq_off+61 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\theta$ - Integer, Parameter :: curl_j_cross_b_phi = curl_mom_eq_off+62 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\phi$ + Integer, Parameter :: curl_j_cross_b_r = curl_mom_eq_off+69 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_r$ + Integer, Parameter :: curl_j_cross_b_theta = curl_mom_eq_off+70 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\theta$ + Integer, Parameter :: curl_j_cross_b_phi = curl_mom_eq_off+71 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\phi$ - Integer, Parameter :: curl_j_cross_b_r_squared = curl_mom_eq_off+63 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_r\right)^2$ - Integer, Parameter :: curl_j_cross_b_theta_squared = curl_mom_eq_off+64 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\theta\right)^2$ - Integer, Parameter :: curl_j_cross_b_phi_squared = curl_mom_eq_off+65 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\phi\right)^2$ + Integer, Parameter :: curl_j_cross_b_r_squared = curl_mom_eq_off+72 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_r\right)^2$ + Integer, Parameter :: curl_j_cross_b_theta_squared = curl_mom_eq_off+73 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\theta\right)^2$ + Integer, Parameter :: curl_j_cross_b_phi_squared = curl_mom_eq_off+74 ! :tex: $\left(\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B)\times\boldsymbol B\right)\right)\right]_\phi\right)^2$ -! Integer, Parameter :: curl_jp_cross_bm_r = curl_mom_eq_off+66 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_r$ -! Integer, Parameter :: curl_jp_cross_bm_theta = curl_mom_eq_off+67 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_\theta$ -! Integer, Parameter :: curl_jp_cross_bm_phi = curl_mom_eq_off+68 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_\phi$ +! Integer, Parameter :: curl_jp_cross_bm_r = curl_mom_eq_off+75 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_r$ +! Integer, Parameter :: curl_jp_cross_bm_theta = curl_mom_eq_off+76 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_\theta$ +! Integer, Parameter :: curl_jp_cross_bm_phi = curl_mom_eq_off+77 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\overline{\boldsymbol B}\right)\right)\right]_\phi$ -! Integer, Parameter :: curl_jm_cross_bp_r = curl_mom_eq_off+69 ! :tex: $\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_r$ -! Integer, Parameter :: curl_jm_cross_bp_theta = curl_mom_eq_off+70 ! :tex: $\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_\theta$ -! Integer, Parameter :: curl_jm_cross_bp_phi = curl_mom_eq_off+71 ! :tex: $\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_\phi$ +! Integer, Parameter :: curl_jm_cross_bp_r = curl_mom_eq_off+78 ! :tex: $\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_r$ +! Integer, Parameter :: curl_jm_cross_bp_theta = curl_mom_eq_off+79 ! :tex: $\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_\theta$ +! Integer, Parameter :: curl_jm_cross_bp_phi = curl_mom_eq_off+80 ! :tex: $\left[\boldsymbol{\nabla}\times c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\boldsymbol B'\right)\right]_\phi$ -! Integer, Parameter :: curl_jm_cross_bm_r = curl_mom_eq_off+72 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_r$ -! Integer, Parameter :: curl_jm_cross_bm_theta = curl_mom_eq_off+73 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_\theta$ -! Integer, Parameter :: curl_jm_cross_bm_phi = curl_mom_eq_off+74 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_\phi$ +! Integer, Parameter :: curl_jm_cross_bm_r = curl_mom_eq_off+81 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_r$ +! Integer, Parameter :: curl_jm_cross_bm_theta = curl_mom_eq_off+82 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_\theta$ +! Integer, Parameter :: curl_jm_cross_bm_phi = curl_mom_eq_off+83 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\boldsymbol B')\times\boldsymbol B'\right)\right)\right]_\phi$ -! Integer, Parameter :: curl_jp_cross_bp_r = curl_mom_eq_off+74 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_r$ -! Integer, Parameter :: curl_jp_cross_bp_theta = curl_mom_eq_off+75 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_\theta$ -! Integer, Parameter :: curl_jp_cross_bp_phi = curl_mom_eq_off+76 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_\phi$ +! Integer, Parameter :: curl_jp_cross_bp_r = curl_mom_eq_off+84 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_r$ +! Integer, Parameter :: curl_jp_cross_bp_theta = curl_mom_eq_off+85 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_\theta$ +! Integer, Parameter :: curl_jp_cross_bp_phi = curl_mom_eq_off+86 ! :tex: $\left[\boldsymbol{\nabla}\times\left(c_4\left((\boldsymbol{\nabla}\times\overline{\boldsymbol B})\times\overline{\boldsymbol B}\right)\right)\right]_\phi$ From 07b8cd3d5a788d94db4b705842e906b27a07d4b7 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Wed, 22 Jul 2026 19:18:03 -0400 Subject: [PATCH 12/13] A logic fix for the squared versions of the curl of the viscous forces. --- src/Diagnostics/Diagnostics_Curl_Momentum.F90 | 6 +++--- src/Diagnostics/Diagnostics_Linear_Forces.F90 | 10 ++++++++-- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/src/Diagnostics/Diagnostics_Curl_Momentum.F90 b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 index 07a4103f2..5d92f069b 100644 --- a/src/Diagnostics/Diagnostics_Curl_Momentum.F90 +++ b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 @@ -630,17 +630,17 @@ Subroutine Initialize_Grad_Viscous_Force() compute_vforce_i_dj = .false. vfoff = 0 - If (sometimes_compute(curl_viscous_force_r)) Then + If (sometimes_compute(curl_viscous_force_r) .or. sometimes_compute(curl_viscous_force_r_squared)) Then compute_vforce_i_dj(2,3) = .true. compute_vforce_i_dj(3,2) = .true. Endif - If (sometimes_compute(curl_viscous_force_theta)) Then + If (sometimes_compute(curl_viscous_force_theta) .or. sometimes_compute(curl_viscous_force_theta_squared)) Then compute_vforce_i_dj(1,3) = .true. compute_vforce_i_dj(3,1) = .true. Endif - If (sometimes_compute(curl_viscous_force_phi)) Then + If (sometimes_compute(curl_viscous_force_phi) .or. sometimes_compute(curl_viscous_force_phi_squared)) Then compute_vforce_i_dj(1,2) = .true. compute_vforce_i_dj(2,1) = .true. Endif diff --git a/src/Diagnostics/Diagnostics_Linear_Forces.F90 b/src/Diagnostics/Diagnostics_Linear_Forces.F90 index e93a5348d..be8ceb48e 100644 --- a/src/Diagnostics/Diagnostics_Linear_Forces.F90 +++ b/src/Diagnostics/Diagnostics_Linear_Forces.F90 @@ -520,7 +520,9 @@ Subroutine Initialize_Viscous_Force() 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 nvf = nvf+1 vf_r = nvf @@ -529,7 +531,9 @@ Subroutine Initialize_Viscous_Force() If (compute_quantity(viscous_force_theta) .or. & compute_quantity(visc_work) .or. & compute_quantity(curl_viscous_force_r) .or. & - compute_quantity(curl_viscous_force_phi)) Then + compute_quantity(curl_viscous_force_r_squared) .or. & + compute_quantity(curl_viscous_force_phi) .or. & + compute_quantity(curl_viscous_force_phi_squared)) Then nvf = nvf + 1 vf_t = nvf Endif @@ -537,7 +541,9 @@ Subroutine Initialize_Viscous_Force() If (compute_quantity(viscous_force_phi) .or. & compute_quantity(visc_work) .or. & compute_quantity(curl_viscous_force_r) .or. & - compute_quantity(curl_viscous_force_theta)) Then + compute_quantity(curl_viscous_force_r_squared) .or. & + compute_quantity(curl_viscous_force_theta) .or. & + compute_quantity(curl_viscous_force_theta_squared)) Then nvf = nvf + 1 vf_p = nvf Endif From aa96b560316e3fe4b6ab36fc7fbe3fd484cf3464 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Fri, 24 Jul 2026 16:49:33 -0400 Subject: [PATCH 13/13] A series of fixes for the curl of momentum diagnostics. The logic was incorrect for when to evaluate the viscous forces as they now need to be evaluated when just the curl is requested. Also fixing the logic for outputting of other curl quantities that previously required both the value and the square to be included. The curl of the momentum terms was never called. Adding an input file that doesn't currently run a test but it is what I've been using to debug these things. --- src/Diagnostics/Diagnostics_Curl_Momentum.F90 | 634 ++++++++---------- src/Diagnostics/Diagnostics_Interface.F90 | 5 +- src/Diagnostics/Diagnostics_Linear_Forces.F90 | 135 ++-- .../Diagnostics_Second_Derivatives.F90 | 58 +- tests/curl_mom_diagnostics/main_input | 71 ++ 5 files changed, 486 insertions(+), 417 deletions(-) create mode 100644 tests/curl_mom_diagnostics/main_input diff --git a/src/Diagnostics/Diagnostics_Curl_Momentum.F90 b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 index 5d92f069b..a3cce30ec 100644 --- a/src/Diagnostics/Diagnostics_Curl_Momentum.F90 +++ b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 @@ -53,35 +53,32 @@ Subroutine Compute_Curl_Advection_Force(buffer) !!!!!!!!!!!!!!!!!!!!!!!!! Advection Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! If (compute_quantity(curl_v_grad_v_r) .or. compute_quantity(curl_v_grad_v_r_squared)) Then - If (compute_quantity(curl_v_grad_v_r)) Then - DO_PSI - qty(PSI) = one_over_r(r) * ref%density(r) * ( & - buffer(PSI,vr) * one_over_r(r) * DDBUFF(PSI,dvpdrdt) + & - buffer(PSI,dvrdt) * one_over_r(r) * buffer(PSI,dvpdr) + & - buffer(PSI,vtheta) * one_over_r(r) * DDBUFF(PSI,dvpdtdt) + & - buffer(PSI,dvtdt) * one_over_r(r) * buffer(PSI,dvpdt) + & - buffer(PSI,vphi) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvpdtdp) - & - csctheta(t) * csctheta(t) * costheta(t) * buffer(PSI,dvpdp) + & - buffer(PSI,dvrdt) + & - cottheta(t) * buffer(PSI,dvtdt) - & - csctheta(t) * csctheta(t) * buffer(PSI,vtheta)) + & - buffer(PSI,dvpdt) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvpdp) + & - buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) + & - one_over_r(r) * costheta(t) * csctheta(t) * ref%density(r) * (buffer(PSI,vr) * & - one_over_r(r) * buffer(PSI,dvpdr) + buffer(PSI,vtheta) * one_over_r(r) * & - buffer(PSI,dvpdt) + buffer(PSI,vphi) * one_over_r(r) * (csctheta(t) * & - buffer(PSI,dvpdp) + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) - & - one_over_r(r) * csctheta(t) * ref%density(r) * (buffer(PSI,vr) * DDBUFF(PSI,dvtdrdp) + & - buffer(PSI,dvrdp) * buffer(PSI,dvtdr) + buffer(PSI,vtheta) * one_over_r(r) * & - (DDBUFF(PSI,dvtdtdp) + buffer(PSI,dvrdp)) + buffer(PSI,dvtdp) * one_over_r(r) * & - (buffer(PSI,dvtdt) + buffer(PSI,vr)) + buffer(PSI,vphi) * one_over_r(r) * & - (csctheta(t) * DDBUFF(PSI,dvtdpdp) - cottheta(t) * buffer(PSI,dvpdp)) + & - buffer(PSI,dvpdp) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvtdp) - & - cottheta(t) * buffer(PSI,vphi))) - END_DO - Call Add_Quantity(qty) - Endif - + DO_PSI + qty(PSI) = one_over_r(r) * ref%density(r) * ( & + buffer(PSI,vr) * one_over_r(r) * DDBUFF(PSI,dvpdrdt) + & + buffer(PSI,dvrdt) * one_over_r(r) * buffer(PSI,dvpdr) + & + buffer(PSI,vtheta) * one_over_r(r) * DDBUFF(PSI,dvpdtdt) + & + buffer(PSI,dvtdt) * one_over_r(r) * buffer(PSI,dvpdt) + & + buffer(PSI,vphi) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvpdtdp) - & + csctheta(t) * csctheta(t) * costheta(t) * buffer(PSI,dvpdp) + & + buffer(PSI,dvrdt) + & + cottheta(t) * buffer(PSI,dvtdt) - & + csctheta(t) * csctheta(t) * buffer(PSI,vtheta)) + & + buffer(PSI,dvpdt) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvpdp) + & + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) + & + one_over_r(r) * costheta(t) * csctheta(t) * ref%density(r) * (buffer(PSI,vr) * & + one_over_r(r) * buffer(PSI,dvpdr) + buffer(PSI,vtheta) * one_over_r(r) * & + buffer(PSI,dvpdt) + buffer(PSI,vphi) * one_over_r(r) * (csctheta(t) * & + buffer(PSI,dvpdp) + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) - & + one_over_r(r) * csctheta(t) * ref%density(r) * (buffer(PSI,vr) * DDBUFF(PSI,dvtdrdp) + & + buffer(PSI,dvrdp) * buffer(PSI,dvtdr) + buffer(PSI,vtheta) * one_over_r(r) * & + (DDBUFF(PSI,dvtdtdp) + buffer(PSI,dvrdp)) + buffer(PSI,dvtdp) * one_over_r(r) * & + (buffer(PSI,dvtdt) + buffer(PSI,vr)) + buffer(PSI,vphi) * one_over_r(r) * & + (csctheta(t) * DDBUFF(PSI,dvtdpdp) - cottheta(t) * buffer(PSI,dvpdp)) + & + buffer(PSI,dvpdp) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvtdp) - & + cottheta(t) * buffer(PSI,vphi))) + END_DO + If (compute_quantity(curl_v_grad_v_r)) Call Add_Quantity(qty) If (compute_quantity(curl_v_grad_v_r_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -92,38 +89,36 @@ Subroutine Compute_Curl_Advection_Force(buffer) Endif If (compute_quantity(curl_v_grad_v_theta) .or. compute_quantity(curl_v_grad_v_theta_squared)) Then - If (compute_quantity(curl_v_grad_v_theta)) Then - DO_PSI - qty(PSI) = one_over_r(r) * csctheta(t) * ref%density(r) * ( buffer(PSI,vr) * DDBUFF(PSI,dvrdrdp) + & - buffer(PSI,dvrdp) * buffer(PSI,dvrdr) + buffer(PSI,vtheta) * one_over_r(r) * & - (DDBUFF(PSI,dvrdtdp) - buffer(PSI,dvtdp)) + buffer(PSI,dvtdp) * one_over_r(r) * & - (buffer(PSI,dvrdt) - buffer(PSI,vtheta)) + buffer(PSI,vphi) * one_over_r(r) * & - (csctheta(t) * DDBUFF(PSI,dvrdpdp) - buffer(PSI,dvpdp)) + buffer(PSI,dvpdp) * & - one_over_r(r) * ( csctheta(t) * buffer(PSI,dvrdp) - buffer(PSI,vphi))) - & - one_over_r(r) * ref%density(r) * (buffer(PSI,vr) * one_over_r(r) * buffer(PSI,dvpdr) + & - buffer(PSI,vtheta) * one_over_r(r) * buffer(PSI,dvpdt) + buffer(PSI,vphi) * one_over_r(r) * & - (csctheta(t) * buffer(PSI,dvpdp) + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) - & - ref%density(r) * (buffer(PSI,vr) * one_over_r(r) * DDBUFF(PSI,dvpdrdr) + & - buffer(PSI,dvrdr) * one_over_r(r) * buffer(PSI,dvpdr) - & - buffer(PSI,vr) * one_over_r(r) * one_over_r(r) * buffer(PSI,dvpdr) + & - buffer(PSI,vtheta) * one_over_r(r) * DDBUFF(PSI,dvpdrdt) + & - buffer(PSI,dvtdr) * one_over_r(r) * buffer(PSI,dvpdt) - & - buffer(PSI,vtheta) * one_over_r(r) * one_over_r(r) * buffer(PSI,dvpdt) + & - buffer(PSI,vtheta) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvpdrdp) + & - buffer(PSI,dvrdr) + & - cottheta(t) * buffer(PSI,dvtdr)) + & - buffer(PSI,dvpdr) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvpdp) + & - buffer(PSI,vr) + & - cottheta(t) * buffer(PSI,vtheta)) - & - buffer(PSI,vphi) * one_over_r(r) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvpdp) + & - buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) - & - ref%dlnrho(r) * & - ref%density(r) * (buffer(PSI,vr) * one_over_r(r) * buffer(PSI,dvpdr) + & - buffer(PSI,vtheta) * one_over_r(r) * buffer(PSI,dvpdt) + buffer(PSI,vphi) * one_over_r(r) * & - (csctheta(t) * buffer(PSI,dvpdp) + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) - END_DO - Call Add_Quantity(qty) - Endif + DO_PSI + qty(PSI) = one_over_r(r) * csctheta(t) * ref%density(r) * ( buffer(PSI,vr) * DDBUFF(PSI,dvrdrdp) + & + buffer(PSI,dvrdp) * buffer(PSI,dvrdr) + buffer(PSI,vtheta) * one_over_r(r) * & + (DDBUFF(PSI,dvrdtdp) - buffer(PSI,dvtdp)) + buffer(PSI,dvtdp) * one_over_r(r) * & + (buffer(PSI,dvrdt) - buffer(PSI,vtheta)) + buffer(PSI,vphi) * one_over_r(r) * & + (csctheta(t) * DDBUFF(PSI,dvrdpdp) - buffer(PSI,dvpdp)) + buffer(PSI,dvpdp) * & + one_over_r(r) * ( csctheta(t) * buffer(PSI,dvrdp) - buffer(PSI,vphi))) - & + one_over_r(r) * ref%density(r) * (buffer(PSI,vr) * one_over_r(r) * buffer(PSI,dvpdr) + & + buffer(PSI,vtheta) * one_over_r(r) * buffer(PSI,dvpdt) + buffer(PSI,vphi) * one_over_r(r) * & + (csctheta(t) * buffer(PSI,dvpdp) + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) - & + ref%density(r) * (buffer(PSI,vr) * one_over_r(r) * DDBUFF(PSI,dvpdrdr) + & + buffer(PSI,dvrdr) * one_over_r(r) * buffer(PSI,dvpdr) - & + buffer(PSI,vr) * one_over_r(r) * one_over_r(r) * buffer(PSI,dvpdr) + & + buffer(PSI,vtheta) * one_over_r(r) * DDBUFF(PSI,dvpdrdt) + & + buffer(PSI,dvtdr) * one_over_r(r) * buffer(PSI,dvpdt) - & + buffer(PSI,vtheta) * one_over_r(r) * one_over_r(r) * buffer(PSI,dvpdt) + & + buffer(PSI,vtheta) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvpdrdp) + & + buffer(PSI,dvrdr) + & + cottheta(t) * buffer(PSI,dvtdr)) + & + buffer(PSI,dvpdr) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvpdp) + & + buffer(PSI,vr) + & + cottheta(t) * buffer(PSI,vtheta)) - & + buffer(PSI,vphi) * one_over_r(r) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvpdp) + & + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) - & + ref%dlnrho(r) * & + ref%density(r) * (buffer(PSI,vr) * one_over_r(r) * buffer(PSI,dvpdr) + & + buffer(PSI,vtheta) * one_over_r(r) * buffer(PSI,dvpdt) + buffer(PSI,vphi) * one_over_r(r) * & + (csctheta(t) * buffer(PSI,dvpdp) + buffer(PSI,vr) + cottheta(t) * buffer(PSI,vtheta))) + END_DO + If (compute_quantity(curl_v_grad_v_theta)) Call Add_Quantity(qty) If (compute_quantity(curl_v_grad_v_theta_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -134,41 +129,39 @@ Subroutine Compute_Curl_Advection_Force(buffer) Endif If (compute_quantity(curl_v_grad_v_phi) .or. compute_quantity(curl_v_grad_v_phi_squared)) Then - If (compute_quantity(curl_v_grad_v_theta)) Then - DO_PSI - qty(PSI) = ref%density(r) * (buffer(PSI,vr) * DDBUFF(PSI,dvtdrdr) + & - buffer(PSI,dvrdr) * buffer(PSI,dvtdr) + & - buffer(PSI,vtheta) * one_over_r(r) * (DDBUFF(PSI,dvtdrdt) + & - buffer(PSI,dvrdr)) + & - buffer(PSI,dvtdr) * one_over_r(r) * (buffer(PSI,dvtdt) + & - buffer(PSI,vr)) - & - buffer(PSI,vtheta) * one_over_r(r) * one_over_r(r) * (buffer(PSI,dvtdt) + & - buffer(PSI,vr)) + & - buffer(PSI,vtheta) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvtdrdp) - & - cottheta(t) * buffer(PSI,dvpdr)) + & - buffer(PSI,dvpdr) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvtdp) - & - cottheta(t) * buffer(PSI,vphi)) - & - buffer(PSI,vphi) * one_over_r(r) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvtdp) - & - cottheta(t) * buffer(PSI,vphi))) + & - ref%dlnrho(r) * ref%density(r) * (buffer(PSI,vr) * buffer(PSI,dvtdr) + buffer(PSI,vtheta) * & - one_over_r(r) * (buffer(PSI,dvtdt) + buffer(PSI,vr)) + buffer(PSI,vphi) * one_over_r(r) * & - (csctheta(t) * buffer(PSI,dvtdp) - cottheta(t) * buffer(PSI,vphi))) + & - one_over_r(r) * ref%density(r) * (buffer(PSI,vr) * buffer(PSI,dvtdr) + buffer(PSI,vtheta) * & - one_over_r(r) * (buffer(PSI,dvtdt) + buffer(PSI,vr)) + buffer(PSI,vphi) * one_over_r(r) * & - (csctheta(t) * buffer(PSI,dvtdp) - cottheta(t) * buffer(PSI,vphi))) - & - one_over_r(r) * ref%density(r) * (buffer(PSI,vr) * DDBUFF(PSI,dvrdrdt) + & - buffer(PSI,dvrdt) * buffer(PSI,dvrdr) + buffer(PSI,vtheta) * one_over_r(r) * & - (DDBUFF(PSI,dvrdtdt) + buffer(PSI,dvtdt)) + & - buffer(PSI,dvtdt) * one_over_r(r) * (buffer(PSI,dvrdt) - buffer(PSI,vtheta)) + & - buffer(PSI,vphi) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvrdtdp) - & - csctheta(t) * csctheta(t) * costheta(t) * buffer(PSI,dvrdp) - & - buffer(PSI,dvpdt)) + & - buffer(PSI,dvpdt) * one_over_r(r) * & - (csctheta(t) * buffer(PSI,dvrdp) - buffer(PSI,vphi))) - - END_DO - Call Add_Quantity(qty) - Endif + DO_PSI + qty(PSI) = ref%density(r) * (buffer(PSI,vr) * DDBUFF(PSI,dvtdrdr) + & + buffer(PSI,dvrdr) * buffer(PSI,dvtdr) + & + buffer(PSI,vtheta) * one_over_r(r) * (DDBUFF(PSI,dvtdrdt) + & + buffer(PSI,dvrdr)) + & + buffer(PSI,dvtdr) * one_over_r(r) * (buffer(PSI,dvtdt) + & + buffer(PSI,vr)) - & + buffer(PSI,vtheta) * one_over_r(r) * one_over_r(r) * (buffer(PSI,dvtdt) + & + buffer(PSI,vr)) + & + buffer(PSI,vtheta) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvtdrdp) - & + cottheta(t) * buffer(PSI,dvpdr)) + & + buffer(PSI,dvpdr) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvtdp) - & + cottheta(t) * buffer(PSI,vphi)) - & + buffer(PSI,vphi) * one_over_r(r) * one_over_r(r) * (csctheta(t) * buffer(PSI,dvtdp) - & + cottheta(t) * buffer(PSI,vphi))) + & + ref%dlnrho(r) * ref%density(r) * (buffer(PSI,vr) * buffer(PSI,dvtdr) + buffer(PSI,vtheta) * & + one_over_r(r) * (buffer(PSI,dvtdt) + buffer(PSI,vr)) + buffer(PSI,vphi) * one_over_r(r) * & + (csctheta(t) * buffer(PSI,dvtdp) - cottheta(t) * buffer(PSI,vphi))) + & + one_over_r(r) * ref%density(r) * (buffer(PSI,vr) * buffer(PSI,dvtdr) + buffer(PSI,vtheta) * & + one_over_r(r) * (buffer(PSI,dvtdt) + buffer(PSI,vr)) + buffer(PSI,vphi) * one_over_r(r) * & + (csctheta(t) * buffer(PSI,dvtdp) - cottheta(t) * buffer(PSI,vphi))) - & + one_over_r(r) * ref%density(r) * (buffer(PSI,vr) * DDBUFF(PSI,dvrdrdt) + & + buffer(PSI,dvrdt) * buffer(PSI,dvrdr) + buffer(PSI,vtheta) * one_over_r(r) * & + (DDBUFF(PSI,dvrdtdt) + buffer(PSI,dvtdt)) + & + buffer(PSI,dvtdt) * one_over_r(r) * (buffer(PSI,dvrdt) - buffer(PSI,vtheta)) + & + buffer(PSI,vphi) * one_over_r(r) * (csctheta(t) * DDBUFF(PSI,dvrdtdp) - & + csctheta(t) * csctheta(t) * costheta(t) * buffer(PSI,dvrdp) - & + buffer(PSI,dvpdt)) + & + buffer(PSI,dvpdt) * one_over_r(r) * & + (csctheta(t) * buffer(PSI,dvrdp) - buffer(PSI,vphi))) + + END_DO + If (compute_quantity(curl_v_grad_v_theta)) Call Add_Quantity(qty) If (compute_quantity(curl_v_grad_v_phi_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -189,14 +182,11 @@ Subroutine Compute_Curl_Buoyancy_Force(buffer) If (compute_quantity(curl_buoyancy_force_theta) .or. compute_quantity(curl_buoyancy_force_theta_squared)) Then - If (compute_quantity(curl_buoyancy_force_theta)) Then - DO_PSI - qty(PSI) = ref%Buoyancy_Coeff(r) * (csctheta(t) * & - radius(r) * buffer(PSI,dtdp)) ! since dtdp = (1/r)*d(temperature or entropy)/dphi - END_DO - Call Add_Quantity(qty) - Endif - + DO_PSI + qty(PSI) = ref%Buoyancy_Coeff(r) * (csctheta(t) * & + radius(r) * buffer(PSI,dtdp)) ! since dtdp = (1/r)*d(temperature or entropy)/dphi + END_DO + If (compute_quantity(curl_buoyancy_force_theta)) Call Add_Quantity(qty) If (compute_quantity(curl_buoyancy_force_theta_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -207,13 +197,10 @@ Subroutine Compute_Curl_Buoyancy_Force(buffer) Endif If (compute_quantity(curl_buoyancy_force_phi) .or. compute_quantity(curl_buoyancy_force_phi_squared)) Then - If (compute_quantity(curl_buoyancy_force_phi)) Then - DO_PSI - qty(PSI) = -ref%Buoyancy_Coeff(r) * ( radius(r) * buffer(PSI,dtdt)) ! since dtdt = (1/r)*d(temperature or $ - END_DO - Call Add_Quantity(qty) - Endif - + DO_PSI + qty(PSI) = -ref%Buoyancy_Coeff(r) * ( radius(r) * buffer(PSI,dtdt)) ! since dtdt = (1/r)*d(temperature or $ + END_DO + If (compute_quantity(curl_buoyancy_force_phi)) Call Add_Quantity(qty) If (compute_quantity(curl_buoyancy_force_phi_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -233,57 +220,54 @@ Subroutine Compute_Curl_Magnetic_Force(buffer) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Magnetic Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! If (compute_quantity(curl_j_cross_b_r) .or. compute_quantity(curl_j_cross_b_r_squared)) Then - If (compute_quantity(curl_j_cross_b_r)) Then - DO_PSI - qty(PSI) = ref%Lorentz_Coeff*(one_over_r(r) * one_over_r(r) * buffer(PSI,dbtdt) * buffer(PSI,dbpdt) + & - one_over_r(r) * one_over_r(r) * buffer(PSI,btheta) * DDBUFF(PSI,dbpdtdt) - & - one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * buffer(PSI,btheta) * & - buffer(PSI,bphi) + & - one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,dbtdt) * & - buffer(PSI,bphi) + & - one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,btheta) * & - buffer(PSI,dbpdt) - & - one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * & - buffer(PSI,dbtdt) * buffer(PSI,dbtdp) - & - one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * & - buffer(PSI,btheta) * DDBUFF(PSI,dbtdtdp) + & - one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * buffer(PSI,btheta) * & - buffer(PSI,dbtdp) + & - one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * & - buffer(PSI,br) * buffer(PSI,dbrdp) - & - one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,dbrdt) * buffer(PSI,dbrdp) - & - one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,br) * DDBUFF(PSI,dbrdtdp) + & - one_over_r(r) * one_over_r(r) * buffer(PSI,dbrdt) * buffer(PSI,bphi) + & - one_over_r(r) * one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbpdt) + & - one_over_r(r) * buffer(PSI,dbrdt) * buffer(PSI,dbpdr) + & - one_over_r(r) * buffer(PSI,br) * DDBUFF(PSI,dbpdrdt) + & !(first part ended) - one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & - buffer(PSI,btheta) * buffer(PSI,dbpdt) + & - one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * costheta(t) * & - buffer(PSI,btheta) * buffer(PSI,bphi) - & - one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & - one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * & - buffer(PSI,btheta) * buffer(PSI,dbtdp) - & - one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * & - buffer(PSI,br) * buffer(PSI,dbrdp) + & - one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & - buffer(PSI,br) * buffer(PSI,bphi) + & - one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,br) * buffer(PSI,dbpdr) - & !(second part ended) - one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbtdp) - & - one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,dbrdp) - & - buffer(PSI,dbrdp) * buffer(PSI,dbtdr) - & - buffer(PSI,br) * DDBUFF(PSI,dbtdrdp) + & - buffer(PSI,dbrdp) * one_over_r(r) * buffer(PSI,dbrdt) + & - buffer(PSI,br) * one_over_r(r) * DDBUFF(PSI,dbrdtdp) + & - one_over_r(r) * buffer(PSI,dbpdp) * buffer(PSI,dbpdt) + & - one_over_r(r) * buffer(PSI,bphi) * DDBUFF(PSI,dbpdtdp) + & - 2 * one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,bphi) * buffer(PSI,dbpdp) - & - buffer(PSI,dbpdp) * buffer(PSI,dbtdp) - & - buffer(PSI,bphi) * DDBUFF(PSI,dbtdpdp)) - END_DO - Call Add_Quantity(qty) - Endif - + DO_PSI + qty(PSI) = ref%Lorentz_Coeff*(one_over_r(r) * one_over_r(r) * buffer(PSI,dbtdt) * buffer(PSI,dbpdt) + & + one_over_r(r) * one_over_r(r) * buffer(PSI,btheta) * DDBUFF(PSI,dbpdtdt) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * buffer(PSI,btheta) * & + buffer(PSI,bphi) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,dbtdt) * & + buffer(PSI,bphi) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,btheta) * & + buffer(PSI,dbpdt) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * & + buffer(PSI,dbtdt) * buffer(PSI,dbtdp) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * & + buffer(PSI,btheta) * DDBUFF(PSI,dbtdtdp) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * buffer(PSI,btheta) * & + buffer(PSI,dbtdp) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * & + buffer(PSI,br) * buffer(PSI,dbrdp) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,dbrdt) * buffer(PSI,dbrdp) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,br) * DDBUFF(PSI,dbrdtdp) + & + one_over_r(r) * one_over_r(r) * buffer(PSI,dbrdt) * buffer(PSI,bphi) + & + one_over_r(r) * one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbpdt) + & + one_over_r(r) * buffer(PSI,dbrdt) * buffer(PSI,dbpdr) + & + one_over_r(r) * buffer(PSI,br) * DDBUFF(PSI,dbpdrdt) + & !(first part ended) + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & + buffer(PSI,btheta) * buffer(PSI,dbpdt) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * costheta(t) * & + buffer(PSI,btheta) * buffer(PSI,bphi) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * & + buffer(PSI,btheta) * buffer(PSI,dbtdp) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * & + buffer(PSI,br) * buffer(PSI,dbrdp) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & + buffer(PSI,br) * buffer(PSI,bphi) + & + one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,br) * buffer(PSI,dbpdr) - & !(second part ended) + one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbtdp) - & + one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,dbrdp) - & + buffer(PSI,dbrdp) * buffer(PSI,dbtdr) - & + buffer(PSI,br) * DDBUFF(PSI,dbtdrdp) + & + buffer(PSI,dbrdp) * one_over_r(r) * buffer(PSI,dbrdt) + & + buffer(PSI,br) * one_over_r(r) * DDBUFF(PSI,dbrdtdp) + & + one_over_r(r) * buffer(PSI,dbpdp) * buffer(PSI,dbpdt) + & + one_over_r(r) * buffer(PSI,bphi) * DDBUFF(PSI,dbpdtdp) + & + 2 * one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,bphi) * buffer(PSI,dbpdp) - & + buffer(PSI,dbpdp) * buffer(PSI,dbtdp) - & + buffer(PSI,bphi) * DDBUFF(PSI,dbtdpdp)) + END_DO + If (compute_quantity(curl_j_cross_b_r)) Call Add_Quantity(qty) If (compute_quantity(curl_j_cross_b_r_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -295,50 +279,47 @@ Subroutine Compute_Curl_Magnetic_Force(buffer) If (compute_quantity(curl_j_cross_b_theta) .or. compute_quantity(curl_j_cross_b_theta_squared)) Then - If (compute_quantity(curl_j_cross_b_theta)) Then - DO_PSI - qty(PSI) = ref%Lorentz_Coeff*(one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * & - buffer(PSI,dbpdp) * buffer(PSI,dbrdp) + & - one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * & - buffer(PSI,bphi) * DDBUFF(PSI,dbrdpdp) - & - 2 * one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,bphi) * buffer(PSI,dbpdp) - & - one_over_r(r) * csctheta(t) * buffer(PSI,dbpdp) * buffer(PSI,dbpdr) - & - one_over_r(r) * csctheta(t) * buffer(PSI,bphi) * DDBUFF(PSI,dbpdrdp) - & - 2 * one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,btheta) * buffer(PSI,dbtdp) - & - one_over_r(r) * csctheta(t) * buffer(PSI,dbtdp) * buffer(PSI,dbtdr) - & - one_over_r(r) * csctheta(t) * buffer(PSI,btheta) * DDBUFF(PSI,dbtdrdp) + & - one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,dbtdp) * buffer(PSI,dbrdt) + & - one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,btheta) * DDBUFF(PSI,dbrdtdp) - & !(first part ended) - one_over_r(r) * one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,dbpdt) - & - one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & - buffer(PSI,btheta) * buffer(PSI,bphi) + & - one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,btheta) * buffer(PSI,dbtdp) + & - one_over_r(r) * one_over_r(r) * csctheta(t) * & - buffer(PSI,br) * buffer(PSI,dbrdp) - & - one_over_r(r) * one_over_r(r) * buffer(PSI,br) * buffer(PSI,bphi) - & - one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbpdr) + & !(second part ended) - one_over_r(r) * one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,dbpdt) - & - one_over_r(r) * buffer(PSI,dbtdr) * buffer(PSI,dbpdt) - & - one_over_r(r) * buffer(PSI,btheta) * DDBUFF(PSI,dbpdrdt) + & - one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & - buffer(PSI,btheta) * buffer(PSI,bphi) - & - one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,dbtdr) * buffer(PSI,bphi) - & - one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,btheta) * buffer(PSI,dbpdr) + & - one_over_r(r) * csctheta(t) * buffer(PSI,dbtdr) * buffer(PSI,dbtdp) - & - buffer(PSI,btheta) * one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,dbtdr) + & - one_over_r(r) * csctheta(t) * buffer(PSI,btheta) * DDBUFF(PSI,dbtdrdp) - & - one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,br) * buffer(PSI,dbrdp) + & - one_over_r(r) * csctheta(t) * buffer(PSI,dbrdr) * buffer(PSI,dbrdp) + & - one_over_r(r) * csctheta(t) * buffer(PSI,br) * DDBUFF(PSI,dbrdrdp) + & - one_over_r(r) * one_over_r(r) * buffer(PSI,br) * buffer(PSI,bphi) - & - one_over_r(r) * buffer(PSI,dbrdr) * buffer(PSI,bphi) - & - one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbpdr) - & - buffer(PSI,dbrdr) * buffer(PSI,dbpdr) - & - buffer(PSI,br) * DDBUFF(PSI,dbpdrdr)) - END_DO - Call Add_Quantity(qty) - Endif - + DO_PSI + qty(PSI) = ref%Lorentz_Coeff*(one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * & + buffer(PSI,dbpdp) * buffer(PSI,dbrdp) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * csctheta(t) * & + buffer(PSI,bphi) * DDBUFF(PSI,dbrdpdp) - & + 2 * one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,bphi) * buffer(PSI,dbpdp) - & + one_over_r(r) * csctheta(t) * buffer(PSI,dbpdp) * buffer(PSI,dbpdr) - & + one_over_r(r) * csctheta(t) * buffer(PSI,bphi) * DDBUFF(PSI,dbpdrdp) - & + 2 * one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,btheta) * buffer(PSI,dbtdp) - & + one_over_r(r) * csctheta(t) * buffer(PSI,dbtdp) * buffer(PSI,dbtdr) - & + one_over_r(r) * csctheta(t) * buffer(PSI,btheta) * DDBUFF(PSI,dbtdrdp) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,dbtdp) * buffer(PSI,dbrdt) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,btheta) * DDBUFF(PSI,dbrdtdp) - & !(first part ended) + one_over_r(r) * one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,dbpdt) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & + buffer(PSI,btheta) * buffer(PSI,bphi) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,btheta) * buffer(PSI,dbtdp) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * & + buffer(PSI,br) * buffer(PSI,dbrdp) - & + one_over_r(r) * one_over_r(r) * buffer(PSI,br) * buffer(PSI,bphi) - & + one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbpdr) + & !(second part ended) + one_over_r(r) * one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,dbpdt) - & + one_over_r(r) * buffer(PSI,dbtdr) * buffer(PSI,dbpdt) - & + one_over_r(r) * buffer(PSI,btheta) * DDBUFF(PSI,dbpdrdt) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & + buffer(PSI,btheta) * buffer(PSI,bphi) - & + one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,dbtdr) * buffer(PSI,bphi) - & + one_over_r(r) * csctheta(t) * costheta(t) * buffer(PSI,btheta) * buffer(PSI,dbpdr) + & + one_over_r(r) * csctheta(t) * buffer(PSI,dbtdr) * buffer(PSI,dbtdp) - & + buffer(PSI,btheta) * one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,dbtdr) + & + one_over_r(r) * csctheta(t) * buffer(PSI,btheta) * DDBUFF(PSI,dbtdrdp) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI,br) * buffer(PSI,dbrdp) + & + one_over_r(r) * csctheta(t) * buffer(PSI,dbrdr) * buffer(PSI,dbrdp) + & + one_over_r(r) * csctheta(t) * buffer(PSI,br) * DDBUFF(PSI,dbrdrdp) + & + one_over_r(r) * one_over_r(r) * buffer(PSI,br) * buffer(PSI,bphi) - & + one_over_r(r) * buffer(PSI,dbrdr) * buffer(PSI,bphi) - & + one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbpdr) - & + buffer(PSI,dbrdr) * buffer(PSI,dbpdr) - & + buffer(PSI,br) * DDBUFF(PSI,dbpdrdr)) + END_DO + If (compute_quantity(curl_j_cross_b_theta)) Call Add_Quantity(qty) If (compute_quantity(curl_j_cross_b_theta_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -349,48 +330,45 @@ Subroutine Compute_Curl_Magnetic_Force(buffer) Endif If (compute_quantity(curl_j_cross_b_phi) .or. compute_quantity(curl_j_cross_b_phi_squared)) Then - If (compute_quantity(curl_j_cross_b_phi)) Then - DO_PSI - qty(PSI) = ref%Lorentz_Coeff*(one_over_r(r) * one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,br) + & - one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbtdr) - & - one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbrdt) - & - one_over_r(r) * one_over_r(r) * buffer(PSI,bphi) * buffer(PSI,dbpdt) - & - one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * (buffer(PSI, bphi))**2 + & - one_over_r(r) * buffer(PSI,bphi) * buffer(PSI,dbtdp) - & !(first part ended) - one_over_r(r) * one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,br) + & - one_over_r(r) * buffer(PSI,dbtdr) * buffer(PSI,br) + & - one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,dbrdr) + & - buffer(PSI,dbrdr) * buffer(PSI,dbtdr) + & - buffer(PSI,br) * DDBUFF(PSI,dbtdrdr) - & - buffer(PSI,dbrdr) * buffer(PSI,dbrdt) - & - one_over_r(r) * buffer(PSI,br) * DDBUFF(PSI,dbrdrdt) + & - one_over_r(r) * one_over_r(r) * buffer(PSI,dbrdt) + & - one_over_r(r) * one_over_r(r) * buffer(PSI,bphi) * buffer(PSI,dbpdt) - & - one_over_r(r) * buffer(PSI,dbpdr) * buffer(PSI,dbpdt) - & - one_over_r(r) * buffer(PSI,bphi) * DDBUFF(PSI,dbpdrdt) + & - one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & - (buffer(PSI, bphi))**2 - & - 2 * one_over_r(r) * csctheta(t) * costheta(t) * & - buffer(PSI, bphi) * buffer(PSI,dbpdr) + & - one_over_r(r) * csctheta(t) * buffer(PSI, dbpdr) * buffer(PSI, dbtdp) - & - buffer(PSI, bphi) * one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI, dbtdp) + & - one_over_r(r) * csctheta(t) * buffer(PSI, bphi) * DDBUFF(PSI,dbtdrdp) + & !(second part ended) - one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * & - buffer(PSI, bphi) * buffer(PSI,dbrdp) - & - one_over_r(r) * csctheta(t) * buffer(PSI,dbpdt) * buffer(PSI,dbrdp) - & - one_over_r(r) * csctheta(t) * buffer(PSI,bphi) * DDBUFF(PSI,dbrdtdp) + & - 2 * one_over_r(r) * buffer(PSI,bphi) * buffer(PSI,dbpdt) + & - buffer(PSI,dbpdt) * buffer(PSI,dbpdr) + & - buffer(PSI,bphi) * DDBUFF(PSI,dbpdrdt) + & - 2 * one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,dbtdt) + & - buffer(PSI,dbtdt) * buffer(PSI,dbtdr) + & - buffer(PSI,btheta) * DDBUFF(PSI,dbtdrdt) - & - one_over_r(r) * buffer(PSI,dbtdt) * buffer(PSI,dbrdt) - & - one_over_r(r) * buffer(PSI,btheta) * DDBUFF(PSI,dbrdtdt)) - END_DO - Call Add_Quantity(qty) - Endif - + DO_PSI + qty(PSI) = ref%Lorentz_Coeff*(one_over_r(r) * one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,br) + & + one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbtdr) - & + one_over_r(r) * buffer(PSI,br) * buffer(PSI,dbrdt) - & + one_over_r(r) * one_over_r(r) * buffer(PSI,bphi) * buffer(PSI,dbpdt) - & + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * (buffer(PSI, bphi))**2 + & + one_over_r(r) * buffer(PSI,bphi) * buffer(PSI,dbtdp) - & !(first part ended) + one_over_r(r) * one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,br) + & + one_over_r(r) * buffer(PSI,dbtdr) * buffer(PSI,br) + & + one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,dbrdr) + & + buffer(PSI,dbrdr) * buffer(PSI,dbtdr) + & + buffer(PSI,br) * DDBUFF(PSI,dbtdrdr) - & + buffer(PSI,dbrdr) * buffer(PSI,dbrdt) - & + one_over_r(r) * buffer(PSI,br) * DDBUFF(PSI,dbrdrdt) + & + one_over_r(r) * one_over_r(r) * buffer(PSI,dbrdt) + & + one_over_r(r) * one_over_r(r) * buffer(PSI,bphi) * buffer(PSI,dbpdt) - & + one_over_r(r) * buffer(PSI,dbpdr) * buffer(PSI,dbpdt) - & + one_over_r(r) * buffer(PSI,bphi) * DDBUFF(PSI,dbpdrdt) + & + one_over_r(r) * one_over_r(r) * csctheta(t) * costheta(t) * & + (buffer(PSI, bphi))**2 - & + 2 * one_over_r(r) * csctheta(t) * costheta(t) * & + buffer(PSI, bphi) * buffer(PSI,dbpdr) + & + one_over_r(r) * csctheta(t) * buffer(PSI, dbpdr) * buffer(PSI, dbtdp) - & + buffer(PSI, bphi) * one_over_r(r) * one_over_r(r) * csctheta(t) * buffer(PSI, dbtdp) + & + one_over_r(r) * csctheta(t) * buffer(PSI, bphi) * DDBUFF(PSI,dbtdrdp) + & !(second part ended) + one_over_r(r) * csctheta(t) * csctheta(t) * costheta(t) * & + buffer(PSI, bphi) * buffer(PSI,dbrdp) - & + one_over_r(r) * csctheta(t) * buffer(PSI,dbpdt) * buffer(PSI,dbrdp) - & + one_over_r(r) * csctheta(t) * buffer(PSI,bphi) * DDBUFF(PSI,dbrdtdp) + & + 2 * one_over_r(r) * buffer(PSI,bphi) * buffer(PSI,dbpdt) + & + buffer(PSI,dbpdt) * buffer(PSI,dbpdr) + & + buffer(PSI,bphi) * DDBUFF(PSI,dbpdrdt) + & + 2 * one_over_r(r) * buffer(PSI,btheta) * buffer(PSI,dbtdt) + & + buffer(PSI,dbtdt) * buffer(PSI,dbtdr) + & + buffer(PSI,btheta) * DDBUFF(PSI,dbtdrdt) - & + one_over_r(r) * buffer(PSI,dbtdt) * buffer(PSI,dbrdt) - & + one_over_r(r) * buffer(PSI,btheta) * DDBUFF(PSI,dbrdtdt)) + END_DO + If (compute_quantity(curl_j_cross_b_phi)) Call Add_Quantity(qty) If (compute_quantity(curl_j_cross_b_phi_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -410,16 +388,13 @@ Subroutine Compute_Curl_Coriolis_Force(buffer) !!!!!!!!!!!!!!!!!!!!!!!!!!!! Coriolis Force !!!!!!!!!!!!!!!!!!!!!!!!!! If (compute_quantity(curl_coriolis_force_r) .or. compute_quantity(curl_coriolis_force_r_squared)) Then - If (compute_quantity(curl_coriolis_force_r)) Then - DO_PSI - qty(PSI) = - ref%Coriolis_Coeff * ref%density(r) * one_over_r(r) * & - (-sintheta(t) * buffer(PSI,vtheta) + cottheta(t) * csctheta(t) * buffer(PSI,vtheta) + & - costheta(t) * buffer(PSI,dvtdt) + 2 * costheta(t) * buffer(PSI,vr) + sintheta(t) * & - buffer(PSI,dvrdt) + cottheta(t) * buffer(PSI,dvpdt)) - END_DO - Call Add_Quantity(qty) - Endif - + DO_PSI + qty(PSI) = - ref%Coriolis_Coeff * ref%density(r) * one_over_r(r) * & + (-sintheta(t) * buffer(PSI,vtheta) + cottheta(t) * csctheta(t) * buffer(PSI,vtheta) + & + costheta(t) * buffer(PSI,dvtdt) + 2 * costheta(t) * buffer(PSI,vr) + sintheta(t) * & + buffer(PSI,dvrdt) + cottheta(t) * buffer(PSI,dvpdt)) + END_DO + If (compute_quantity(curl_coriolis_force_r)) Call Add_Quantity(qty) If (compute_quantity(curl_coriolis_force_r_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -430,16 +405,13 @@ Subroutine Compute_Curl_Coriolis_Force(buffer) Endif If (compute_quantity(curl_coriolis_force_theta) .or. compute_quantity(curl_coriolis_force_theta_squared)) Then - If (compute_quantity(curl_coriolis_force_theta)) Then - DO_PSI - qty(PSI) = ref%Coriolis_Coeff * ref%density(r) * (one_over_r(r) * (buffer(PSI,dvpdp) + & - costheta(t) * buffer(PSI,vtheta) + sintheta(t) * buffer(PSI,vr)) + costheta(t) * & - buffer(PSI,dvtdr) + sintheta(t) * buffer(PSI,dvrdr) + ref%dlnrho(r) * costheta(t) * & - buffer(PSI,vtheta) + ref%dlnrho(r) * sintheta(t) * buffer(PSI,vr)) - END_DO - Call Add_Quantity(qty) - Endif - + DO_PSI + qty(PSI) = ref%Coriolis_Coeff * ref%density(r) * (one_over_r(r) * (buffer(PSI,dvpdp) + & + costheta(t) * buffer(PSI,vtheta) + sintheta(t) * buffer(PSI,vr)) + costheta(t) * & + buffer(PSI,dvtdr) + sintheta(t) * buffer(PSI,dvrdr) + ref%dlnrho(r) * costheta(t) * & + buffer(PSI,vtheta) + ref%dlnrho(r) * sintheta(t) * buffer(PSI,vr)) + END_DO + If (compute_quantity(curl_coriolis_force_theta)) Call Add_Quantity(qty) If (compute_quantity(curl_coriolis_force_theta_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -450,15 +422,12 @@ Subroutine Compute_Curl_Coriolis_Force(buffer) Endif If (compute_quantity(curl_coriolis_force_phi) .or. compute_quantity(curl_coriolis_force_phi)) Then - If (compute_quantity(curl_coriolis_force_phi)) Then - DO_PSI - qty(PSI) = ref%Coriolis_Coeff * ref%density(r) * (ref%dlnrho(r) * costheta(t) * buffer(PSI,vphi) + & - costheta(t) * buffer(PSI,dvpdr) - one_over_r(r) * sintheta(t) * buffer(PSI,dvpdr)) - - END_DO - Call Add_Quantity(qty) - Endif + DO_PSI + qty(PSI) = ref%Coriolis_Coeff * ref%density(r) * (ref%dlnrho(r) * costheta(t) * buffer(PSI,vphi) + & + costheta(t) * buffer(PSI,dvpdr) - one_over_r(r) * sintheta(t) * buffer(PSI,dvpdr)) + END_DO + If (compute_quantity(curl_coriolis_force_phi)) Call Add_Quantity(qty) If (compute_quantity(curl_coriolis_force_phi_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -482,14 +451,11 @@ Subroutine Compute_Curl_Pressure_Force(buffer) !!!!!!!!!!!!!!!!!!!!!!!!!!! Pressure Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!! If (compute_quantity(curl_pressure_force_r) .or. compute_quantity(curl_pressure_force_r_squared)) Then - If (compute_quantity(curl_pressure_force_r)) Then - DO_PSI - qty(PSI) = pfactor(r) * & - OneOverRSquared(r) * csctheta(t) * (DDBUFF(PSI,dvpdtdp)-DDBUFF(PSI,dvpdtdp)) - END_DO - Call Add_Quantity(qty) - Endif - + DO_PSI + qty(PSI) = pfactor(r) * & + OneOverRSquared(r) * csctheta(t) * (DDBUFF(PSI,dvpdtdp)-DDBUFF(PSI,dvpdtdp)) + END_DO + If (compute_quantity(curl_pressure_force_r)) Call Add_Quantity(qty) If (compute_quantity(curl_pressure_force_r_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -499,15 +465,12 @@ Subroutine Compute_Curl_Pressure_Force(buffer) Endif If (compute_quantity(curl_pressure_force_theta) .or. compute_quantity(curl_pressure_force_theta_squared)) Then - If (compute_quantity(curl_pressure_force_theta)) Then - DO_PSI - qty(PSI) = pfactor(r) * & - OneOverRSquared(r) * csctheta(t) * (DDBUFF(PSI,dvpdrdp) + & - ref%dlnrho(r) * buffer(PSI,dpdp) + DDBUFF(PSI,dvpdrdp)) - END_DO - Call Add_Quantity(qty) - Endif - + DO_PSI + qty(PSI) = pfactor(r) * & + OneOverRSquared(r) * csctheta(t) * (DDBUFF(PSI,dvpdrdp) + & + ref%dlnrho(r) * buffer(PSI,dpdp) + DDBUFF(PSI,dvpdrdp)) + END_DO + If (compute_quantity(curl_pressure_force_theta)) Call Add_Quantity(qty) If (compute_quantity(curl_pressure_force_theta_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -519,15 +482,12 @@ Subroutine Compute_Curl_Pressure_Force(buffer) If (compute_quantity(curl_pressure_force_phi) .or. compute_quantity(curl_pressure_force_phi_squared)) Then - If (compute_quantity(curl_pressure_force_phi)) Then - DO_PSI - qty(PSI) = pfactor(r) * & - one_over_r(r) * csctheta(t) * (-DDBUFF(PSI,dvpdrdt) + DDBUFF(PSI,dvpdrdt) - & - ref%dlnrho(r) * buffer(PSI,dpdt)) - END_DO - Call Add_Quantity(qty) - Endif - + DO_PSI + qty(PSI) = pfactor(r) * & + one_over_r(r) * csctheta(t) * (-DDBUFF(PSI,dvpdrdt) + DDBUFF(PSI,dvpdrdt) - & + ref%dlnrho(r) * buffer(PSI,dpdt)) + END_DO + If (compute_quantity(curl_pressure_force_phi)) Call Add_Quantity(qty) If (compute_quantity(curl_pressure_force_phi_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -546,14 +506,12 @@ Subroutine Compute_Curl_Viscous_Force(buffer) !!!!!!!!!!!!!!!!!!!!!!!!!!! Viscous Force !!!!!!!!!!!!!!!!!!!!!!!!!!!!! If (compute_quantity(curl_viscous_force_r) .or. compute_quantity(curl_viscous_force_r_squared)) Then - If (compute_quantity(curl_viscous_force_r)) Then - DO_PSI - qty(PSI) = One_Over_R(r)*(VFDBUFF(PSI,dvf_p_dt) + & - cottheta(t)*vforce_buffer(PSI,vfp_r) - & - csctheta(t)*VFDBUFF(PSI,dvf_t_dp)) - END_DO - Call Add_Quantity(qty) - Endif + DO_PSI + qty(PSI) = One_Over_R(r)*(VFDBUFF(PSI,dvf_p_dt) + & + cottheta(t)*vforce_buffer(PSI,vf_p) - & + csctheta(t)*VFDBUFF(PSI,dvf_t_dp)) + END_DO + If (compute_quantity(curl_viscous_force_r)) Call Add_Quantity(qty) If (compute_quantity(curl_viscous_force_r_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -563,14 +521,12 @@ Subroutine Compute_Curl_Viscous_Force(buffer) Endif If (compute_quantity(curl_viscous_force_theta) .or. compute_quantity(curl_viscous_force_theta_squared)) Then - If (compute_quantity(curl_viscous_force_theta)) Then - DO_PSI - qty(PSI) = One_Over_R(r)*(csctheta(t)*VFDBUFF(PSI,dvf_r_dp) - & - vforce_buffer(PSI,vfp_p)) - & - VFDBUFF(PSI,dvf_p_dr) - END_DO - Call Add_Quantity(qty) - Endif + DO_PSI + qty(PSI) = One_Over_R(r)*(csctheta(t)*VFDBUFF(PSI,dvf_r_dp) - & + vforce_buffer(PSI,vf_p)) - & + VFDBUFF(PSI,dvf_p_dr) + END_DO + If (compute_quantity(curl_viscous_force_theta)) Call Add_Quantity(qty) If (compute_quantity(curl_viscous_force_theta_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -579,15 +535,13 @@ Subroutine Compute_Curl_Viscous_Force(buffer) Endif Endif - If (compute_quantity(curl_viscous_force_phi) .or. compute_quantity(curl_viscous_force_phi_squared)) Then - If (compute_quantity(curl_viscous_force_phi)) Then - DO_PSI - qty(PSI) = VFDBUFF(PSI,dvf_t_dr) + & - One_Over_R(r)*(vforce_buffer(PSI,vfp_t) - & - VFDBUFF(PSI,dvf_r_dt)) - END_DO - Call Add_Quantity(qty) - Endif + If (compute_quantity(curl_viscous_force_phi) .or. compute_quantity(curl_viscous_force_phi_squared)) Then + DO_PSI + qty(PSI) = VFDBUFF(PSI,dvf_t_dr) + & + One_Over_R(r)*(vforce_buffer(PSI,vf_t) - & + VFDBUFF(PSI,dvf_r_dt)) + END_DO + If (compute_quantity(curl_viscous_force_phi)) Call Add_Quantity(qty) If (compute_quantity(curl_viscous_force_phi_squared)) Then DO_PSI qty(PSI) = qty(PSI)*qty(PSI) @@ -630,17 +584,20 @@ Subroutine Initialize_Grad_Viscous_Force() compute_vforce_i_dj = .false. vfoff = 0 - If (sometimes_compute(curl_viscous_force_r) .or. sometimes_compute(curl_viscous_force_r_squared)) Then + If (sometimes_compute(curl_viscous_force_r) .or. & + sometimes_compute(curl_viscous_force_r_squared)) Then compute_vforce_i_dj(2,3) = .true. compute_vforce_i_dj(3,2) = .true. Endif - If (sometimes_compute(curl_viscous_force_theta) .or. sometimes_compute(curl_viscous_force_theta_squared)) Then + If (sometimes_compute(curl_viscous_force_theta) .or. & + sometimes_compute(curl_viscous_force_theta_squared)) Then compute_vforce_i_dj(1,3) = .true. compute_vforce_i_dj(3,1) = .true. Endif - If (sometimes_compute(curl_viscous_force_phi) .or. sometimes_compute(curl_viscous_force_phi_squared)) Then + If (sometimes_compute(curl_viscous_force_phi) .or. & + sometimes_compute(curl_viscous_force_phi_squared)) Then compute_vforce_i_dj(1,2) = .true. compute_vforce_i_dj(2,1) = .true. Endif @@ -752,7 +709,7 @@ Subroutine Grad_Viscous_Force() ! load the fields we want to take derivatives of do i = 1, nvffields - d_vforce_buffer%p3b(:,:,:,i) = vforce_buffer(:,:,:,vfdindmap(i, 1)) + d_vforce_buffer%p3b(1:n_phi,:,:,i) = vforce_buffer(:,:,:,vfdindmap(i, 1)) enddo ! transform to Fourier m space @@ -761,15 +718,13 @@ Subroutine Grad_Viscous_Force() ! reform to hybrid rlm space call d_vforce_buffer%reform() ! move to p2b - ! deallocate p3b - call d_vforce_buffer%deconstruct('p3b') - ! allocate spectral buffer and transform call d_vforce_buffer%construct('s2b') call Legendre_Transform(d_vforce_buffer%p2b, d_vforce_buffer%s2b) ! deallocate p2b call d_vforce_buffer%deconstruct('p2b') + d_vforce_buffer%config = 's2b' ! reform call d_vforce_buffer%reform() ! move to p1b @@ -784,6 +739,9 @@ Subroutine Grad_Viscous_Force() d_vforce_buffer%p1a = d_vforce_buffer%p1b end if + d_vforce_buffer%p1b = 0.0 + d_vforce_buffer%config = 'p1a' + ! take d_by_dr ! (and transform back to grid space if in Chebyshev) if (chebyshev) then @@ -802,8 +760,7 @@ Subroutine Grad_Viscous_Force() enddo end if - ! moving back, but first we need grid space in p1a - d_vforce_buffer%config = 'p1a' + ! moving back call d_vforce_buffer%deconstruct('p1b') ! reform and start moving back @@ -819,9 +776,10 @@ Subroutine Grad_Viscous_Force() call d_vforce_buffer%construct('p2a') call Legendre_Transform(d_vforce_buffer%s2a, d_vforce_buffer%p2a) call d_vforce_buffer%deconstruct('s2a') + d_vforce_buffer%config = 'p2a' ! reform - call d_vforce_buffer%reform() ! move to p2a + call d_vforce_buffer%reform() ! move to p3a ! take d_by_dphi derivatives do i = 1, nvffields @@ -837,7 +795,7 @@ Subroutine Grad_Viscous_Force() do i = 1, nvffields if (vfdindmap(i,3) .gt. 0) then DO_PSI - d_vforce_buffer%p3a(PSI,i) = d_vforce_buffer%p3a(PSI,i)*csctheta(t) + d_vforce_buffer%p3a(PSI,vfdindmap(i,3)) = d_vforce_buffer%p3a(PSI,vfdindmap(i,3))*csctheta(t) END_DO end if enddo diff --git a/src/Diagnostics/Diagnostics_Interface.F90 b/src/Diagnostics/Diagnostics_Interface.F90 index 1fe54e1c0..f2e2f5bbb 100755 --- a/src/Diagnostics/Diagnostics_Interface.F90 +++ b/src/Diagnostics/Diagnostics_Interface.F90 @@ -175,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() @@ -225,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) @@ -262,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 @@ -343,6 +342,8 @@ Subroutine Initialize_Diagnostics() Call Initialize_Second_Derivatives() + Call Initialize_Viscous_Force() + Call Initialize_Grad_Viscous_Force() Call Initialize_Diagnostics_Buffer() diff --git a/src/Diagnostics/Diagnostics_Linear_Forces.F90 b/src/Diagnostics/Diagnostics_Linear_Forces.F90 index be8ceb48e..3a443dcc8 100644 --- a/src/Diagnostics/Diagnostics_Linear_Forces.F90 +++ b/src/Diagnostics/Diagnostics_Linear_Forces.F90 @@ -517,82 +517,82 @@ Subroutine Initialize_Viscous_Force() Implicit None nvf = 0 - 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 + 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) .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 + 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) .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 + 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) .or. & - compute_quantity(curl_viscous_pforce_theta) .or. & - compute_quantity(curl_viscous_pforce_phi)) 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) .or. & - compute_quantity(curl_viscous_pforce_r) .or. & - compute_quantity(curl_viscous_pforce_phi)) 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).or. & - compute_quantity(curl_viscous_pforce_r) .or. & - compute_quantity(curl_viscous_pforce_theta)) 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) .or. & - compute_quantity(curl_viscous_mforce_theta) .or. & - compute_quantity(curl_viscous_mforce_phi)) 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).or. & - compute_quantity(curl_viscous_mforce_r) .or. & - compute_quantity(curl_viscous_mforce_phi)) 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).or. & - compute_quantity(curl_viscous_mforce_r) .or. & - compute_quantity(curl_viscous_mforce_theta)) 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 @@ -622,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 @@ -653,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 @@ -679,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) @@ -708,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 @@ -739,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 @@ -765,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) @@ -794,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 @@ -823,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 @@ -850,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) @@ -880,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 diff --git a/src/Diagnostics/Diagnostics_Second_Derivatives.F90 b/src/Diagnostics/Diagnostics_Second_Derivatives.F90 index 27ae6b8ec..562ea4135 100644 --- a/src/Diagnostics/Diagnostics_Second_Derivatives.F90 +++ b/src/Diagnostics/Diagnostics_Second_Derivatives.F90 @@ -85,22 +85,49 @@ Subroutine Init_Derivative_Logic() !////////////////////////////////////////////////////////////////// ! Terms related to viscosity - If (sometimes_compute(viscous_force_r)) compute_vr_dd = .true. - If (sometimes_compute(viscous_pforce_r)) compute_vr_dd = .true. - If (sometimes_compute(viscous_mforce_r)) compute_vr_dd = .true. - - If (sometimes_compute(viscous_force_theta)) compute_vt_dd = .true. - If (sometimes_compute(viscous_pforce_theta)) compute_vt_dd = .true. - If (sometimes_compute(viscous_mforce_theta)) compute_vt_dd = .true. + If (sometimes_compute(visc_work) .or. & + sometimes_compute(viscous_force_r) .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_pforce_r) .or. & + sometimes_compute(curl_viscous_pforce_theta) .or. & + sometimes_compute(curl_viscous_pforce_phi) .or. & + sometimes_compute(viscous_mforce_r) .or. & + sometimes_compute(curl_viscous_mforce_theta) .or. & + sometimes_compute(curl_viscous_mforce_phi)) Then + compute_vr_dd = .true. + Endif - If (sometimes_compute(viscous_force_phi)) compute_vp_dd = .true. - If (sometimes_compute(viscous_pforce_phi)) compute_vp_dd = .true. - If (sometimes_compute(viscous_mforce_phi)) compute_vp_dd = .true. - If (sometimes_compute(visc_work) .or. sometimes_compute(visc_work_pp) & - .or. sometimes_compute(visc_work_mm) ) Then - compute_vr_dd = .true. + If (sometimes_compute(visc_work_pp) .or. & + sometimes_compute(viscous_force_theta) .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) .or. & + sometimes_compute(viscous_pforce_theta) .or. & + sometimes_compute(curl_viscous_pforce_r) .or. & + sometimes_compute(curl_viscous_pforce_phi) .or. & + sometimes_compute(viscous_mforce_theta) .or. & + sometimes_compute(curl_viscous_mforce_r) .or. & + sometimes_compute(curl_viscous_mforce_phi)) Then compute_vt_dd = .true. + Endif + + If (sometimes_compute(visc_work_mm) .or. & + sometimes_compute(viscous_force_phi) .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) .or. & + sometimes_compute(viscous_pforce_phi) .or. & + sometimes_compute(curl_viscous_pforce_r) .or. & + sometimes_compute(curl_viscous_pforce_theta) .or. & + sometimes_compute(viscous_mforce_phi) .or. & + sometimes_compute(curl_viscous_mforce_r) .or. & + sometimes_compute(curl_viscous_mforce_theta)) Then compute_vp_dd = .true. Endif @@ -320,11 +347,6 @@ Subroutine Compute_Second_Derivatives(inbuffer) Real*8, Intent(InOut) :: inbuffer(1:,my_r%min:,my_theta%min:,1:) Type(rmcontainer3D), Allocatable :: ddtemp(:) - - - - - ! Here were compute all second derivatives for N variables ! Outline of the process: ! 1) Intialize p3b space of d2buffer diff --git a/tests/curl_mom_diagnostics/main_input b/tests/curl_mom_diagnostics/main_input new file mode 100644 index 000000000..86c28a307 --- /dev/null +++ b/tests/curl_mom_diagnostics/main_input @@ -0,0 +1,71 @@ +&problemsize_namelist + n_r = 32 + n_theta = 32 + nprow = 2 + npcol = 2 + aspect_ratio = 0.35d0 + shell_depth = 1.0d0 +/ +&numerical_controls_namelist +/ +&physical_controls_namelist + benchmark_mode = 0 + rotation = .True. + magnetism = .true. + viscous_heating = .false. + ohmic_heating = .false. + advect_reference_state = .false. +/ +&temporal_controls_namelist + max_time_step = 1.0d-4 + max_iterations = 100 + checkpoint_interval = 10000 + cflmin = 0.4d0 + cflmax = 0.6d0 +/ +&io_controls_namelist +/ +&output_namelist +globalavg_values = 1339, 1340, 1341 +globalavg_frequency = 10 +globalavg_nrec = 1 + +point_probe_values = 1339, 1340, 1341 +point_probe_r_nrm = 0.0, -1.0 +point_probe_theta_nrm = 0.0, -1.0 +point_probe_phi_nrm = 0.20, -0.30, 0.7, -0.8 +point_probe_frequency = 100 +point_probe_nrec = 1 +/ + +&Boundary_Conditions_Namelist +no_slip_boundaries = .true. +strict_L_Conservation = .false. +dtdr_bottom = 0.0d0 +T_Top = 0.0d0 +T_Bottom = 1.0d0 +fix_tvar_top = .true. +fix_tvar_bottom = .true. +fix_dtdr_bottom = .false. +/ +&Initial_Conditions_Namelist +init_type = 1 +magnetic_init_type = 1 +mag_amp = 1.0d0 +temp_amp = 1.0d1 +temp_w = 0.01d4 +restart_iter = -1 +/ +&Test_Namelist +/ +&Reference_Namelist +Ekman_Number = 1.0d-3 +Rayleigh_Number = 1.0d5 +Prandtl_Number = 1.0d0 +Magnetic_Prandtl_Number = 5.0d0 +reference_type = 1 +heating_type = 0 +gravity_power = 1.0d0 +/ +&Transport_Namelist +/