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..64c4c53e3 --- /dev/null +++ b/doc/source/diagnostic_codes/curl_momentum_equation.rst @@ -0,0 +1,85 @@ +Curl Momentum Equation +==================================================================== + +===================================================================================================================================================================== ====== ======================================= + :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 + :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(\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/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/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 new file mode 100644 index 000000000..a3cce30ec --- /dev/null +++ b/src/Diagnostics/Diagnostics_Curl_Momentum.F90 @@ -0,0 +1,805 @@ +! +! 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 +! . +! + +#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) + 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) + Call Compute_Curl_Viscous_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) + 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 + 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) + 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 + 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) + 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 + 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) + END_DO + Call Add_Quantity(qty) + 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + If (compute_quantity(curl_buoyancy_force_theta) .or. compute_quantity(curl_buoyancy_force_theta_squared)) 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 + 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) + END_DO + Call Add_Quantity(qty) + Endif + + Endif + + If (compute_quantity(curl_buoyancy_force_phi) .or. compute_quantity(curl_buoyancy_force_phi_squared)) Then + 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) + END_DO + Call Add_Quantity(qty) + Endif + + Endif + + 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 + 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) + 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 + 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) + 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 + 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) + END_DO + Call Add_Quantity(qty) + Endif + + Endif + + 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 + 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) + END_DO + Call Add_Quantity(qty) + Endif + + Endif + + If (compute_quantity(curl_coriolis_force_theta) .or. compute_quantity(curl_coriolis_force_theta_squared)) 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 + 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) + END_DO + Call Add_Quantity(qty) + Endif + + Endif + + If (compute_quantity(curl_coriolis_force_phi) .or. 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 + 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) + END_DO + Call Add_Quantity(qty) + Endif + + 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 + 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) + END_DO + Call Add_Quantity(qty) + Endif + Endif + + If (compute_quantity(curl_pressure_force_theta) .or. compute_quantity(curl_pressure_force_theta_squared)) 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 + 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) + END_DO + Call Add_Quantity(qty) + Endif + Endif + + + + If (compute_quantity(curl_pressure_force_phi) .or. compute_quantity(curl_pressure_force_phi_squared)) 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 + 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) + END_DO + Call Add_Quantity(qty) + Endif + Endif + + 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) .or. compute_quantity(curl_viscous_force_r_squared)) Then + 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) + END_DO + Call Add_Quantity(qty) + Endif + Endif + + If (compute_quantity(curl_viscous_force_theta) .or. compute_quantity(curl_viscous_force_theta_squared)) Then + 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) + END_DO + Call Add_Quantity(qty) + Endif + 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) + END_DO + Call Add_Quantity(qty) + Endif + Endif + + End Subroutine Compute_Curl_Viscous_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) .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 + 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 + 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(1:n_phi,:,:,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 + + ! 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 + + ! 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 + + 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 + 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) + 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 + 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') + d_vforce_buffer%config = 'p2a' + + ! reform + call d_vforce_buffer%reform() ! move to p3a + + ! 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,vfdindmap(i,3)) = d_vforce_buffer%p3a(PSI,vfdindmap(i,3))*csctheta(t) + END_DO + end if + enddo + + End Subroutine Grad_Viscous_Force + +End Module Diagnostics_Curl_Momentum diff --git a/src/Diagnostics/Diagnostics_Interface.F90 b/src/Diagnostics/Diagnostics_Interface.F90 index 86469b137..f2e2f5bbb 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 @@ -174,7 +175,6 @@ Subroutine PS_Output(buffer,iteration, current_time) Call ComputeM0(buffer,m0_values) Call Compute_Fluctuations(buffer) - Call Initialize_Viscous_Force() Call Initialize_Mean_Correction() @@ -189,6 +189,10 @@ Subroutine PS_Output(buffer,iteration, current_time) over_n_phi = 1.0d0/dble(n_phi) Call Viscous_Force(buffer) ! Pre-calculate the viscous forces and place them in the vforce_buffer + + if (need_vforce_derivatives) then + Call Grad_Viscous_Force() + endif Call Mean_Correction(buffer) ! Remove ell=0 component from radial and theta forces @@ -220,6 +224,7 @@ Subroutine PS_Output(buffer,iteration, current_time) Call Compute_Angular_Momentum_Balance(buffer) Call Compute_Inertial_Terms(buffer) Call Compute_Linear_Forces(buffer) + Call Compute_Curl_Momentum_Forces(buffer) Call Compute_KE_Flux(buffer) @@ -257,7 +262,6 @@ Subroutine PS_Output(buffer,iteration, current_time) Call d2buffer%deconstruct('p3a') DeAllocate(d2_ell0,d2_m0,d2_fbuffer) ENDIF - Call Finalize_Viscous_Force() Call Finalize_Mean_Correction() Endif ! time_to_output(iteration) End Subroutine PS_Output @@ -338,6 +342,10 @@ Subroutine Initialize_Diagnostics() Call Initialize_Second_Derivatives() + Call Initialize_Viscous_Force() + + Call Initialize_Grad_Viscous_Force() + Call Initialize_Diagnostics_Buffer() End Subroutine Initialize_Diagnostics diff --git a/src/Diagnostics/Diagnostics_Linear_Forces.F90 b/src/Diagnostics/Diagnostics_Linear_Forces.F90 index e4598baed..3a443dcc8 100644 --- a/src/Diagnostics/Diagnostics_Linear_Forces.F90 +++ b/src/Diagnostics/Diagnostics_Linear_Forces.F90 @@ -517,58 +517,82 @@ Subroutine Initialize_Viscous_Force() Implicit None nvf = 0 - If (compute_quantity(viscous_force_r) .or. & - compute_quantity(visc_work) .or. & - compute_quantity(viscous_mforce_r)) Then + If (sometimes_compute(viscous_force_r) .or. & + sometimes_compute(visc_work) .or. & + sometimes_compute(curl_viscous_force_theta) .or. & + sometimes_compute(curl_viscous_force_theta_squared) .or. & + sometimes_compute(curl_viscous_force_phi) .or. & + sometimes_compute(curl_viscous_force_phi_squared) .or. & + sometimes_compute(viscous_mforce_r)) Then nvf = nvf+1 vf_r = nvf Endif - If (compute_quantity(viscous_force_theta) .or. & - compute_quantity(visc_work)) Then + If (sometimes_compute(viscous_force_theta) .or. & + sometimes_compute(visc_work) .or. & + sometimes_compute(curl_viscous_force_r) .or. & + sometimes_compute(curl_viscous_force_r_squared) .or. & + sometimes_compute(curl_viscous_force_phi) .or. & + sometimes_compute(curl_viscous_force_phi_squared)) Then nvf = nvf + 1 vf_t = nvf Endif - If (compute_quantity(viscous_force_phi) .or. & - compute_quantity(visc_work)) Then + If (sometimes_compute(viscous_force_phi) .or. & + sometimes_compute(visc_work) .or. & + sometimes_compute(curl_viscous_force_r) .or. & + sometimes_compute(curl_viscous_force_r_squared) .or. & + sometimes_compute(curl_viscous_force_theta) .or. & + sometimes_compute(curl_viscous_force_theta_squared)) Then nvf = nvf + 1 vf_p = nvf Endif - If (compute_quantity(viscous_pforce_r) .or. & - compute_quantity(visc_work_pp)) Then + If (sometimes_compute(viscous_pforce_r) .or. & + sometimes_compute(visc_work_pp) .or. & + sometimes_compute(curl_viscous_pforce_theta) .or. & + sometimes_compute(curl_viscous_pforce_phi)) Then nvf = nvf+1 vfp_r = nvf Endif - If (compute_quantity(viscous_pforce_theta) .or. & - compute_quantity(visc_work_pp)) Then + If (sometimes_compute(viscous_pforce_theta) .or. & + sometimes_compute(visc_work_pp) .or. & + sometimes_compute(curl_viscous_pforce_r) .or. & + sometimes_compute(curl_viscous_pforce_phi)) Then nvf = nvf + 1 vfp_t = nvf Endif - If (compute_quantity(viscous_pforce_phi) .or. & - compute_quantity(visc_work_pp)) Then + If (sometimes_compute(viscous_pforce_phi) .or. & + sometimes_compute(visc_work_pp) .or. & + sometimes_compute(curl_viscous_pforce_r) .or. & + sometimes_compute(curl_viscous_pforce_theta)) Then nvf = nvf + 1 vfp_p = nvf Endif - If (compute_quantity(viscous_mforce_r) .or. & - compute_quantity(visc_work_mm)) Then + If (sometimes_compute(viscous_mforce_r) .or. & + sometimes_compute(visc_work_mm) .or. & + sometimes_compute(curl_viscous_mforce_theta) .or. & + sometimes_compute(curl_viscous_mforce_phi)) Then nvf = nvf+1 vfm_r = nvf Endif - If (compute_quantity(viscous_mforce_theta) .or. & - compute_quantity(visc_work_mm)) Then + If (sometimes_compute(viscous_mforce_theta) .or. & + sometimes_compute(visc_work_mm) .or. & + sometimes_compute(curl_viscous_mforce_r) .or. & + sometimes_compute(curl_viscous_mforce_phi)) Then nvf = nvf + 1 vfm_t = nvf Endif - If (compute_quantity(viscous_mforce_phi) .or. & - compute_quantity(samom_diffusion) .or. & - compute_quantity(visc_work_mm)) Then + If (sometimes_compute(viscous_mforce_phi) .or. & + sometimes_compute(samom_diffusion) .or. & + sometimes_compute(visc_work_mm) .or. & + sometimes_compute(curl_viscous_mforce_r) .or. & + sometimes_compute(curl_viscous_mforce_theta)) Then nvf = nvf + 1 vfm_p = nvf Endif @@ -598,6 +622,10 @@ Subroutine Viscous_Force(buffer) If (compute_quantity(viscous_force_r) .or. & compute_quantity(visc_work) .or. & + compute_quantity(curl_viscous_force_theta) .or. & + compute_quantity(curl_viscous_force_theta_squared) .or. & + compute_quantity(curl_viscous_force_phi) .or. & + compute_quantity(curl_viscous_force_phi_squared) .or. & compute_quantity(viscous_mforce_r) ) Then DO_PSI @@ -629,7 +657,11 @@ Subroutine Viscous_Force(buffer) !Theta-direction; Full If (compute_quantity(viscous_force_theta) .or. & - compute_quantity(visc_work)) Then + compute_quantity(visc_work) .or. & + compute_quantity(curl_viscous_force_r) .or. & + compute_quantity(curl_viscous_force_r_squared) .or. & + compute_quantity(curl_viscous_force_phi) .or. & + compute_quantity(curl_viscous_force_phi_squared)) Then DO_PSI ! first, compute all the terms multiplied by mu @@ -655,7 +687,11 @@ Subroutine Viscous_Force(buffer) !Phi-direction If (compute_quantity(viscous_force_phi) .or. & - compute_quantity(visc_work)) Then + compute_quantity(visc_work) .or. & + compute_quantity(curl_viscous_force_r) .or. & + compute_quantity(curl_viscous_force_r_squared) .or. & + compute_quantity(curl_viscous_force_theta) .or. & + compute_quantity(curl_viscous_force_theta_squared)) Then DO_PSI del2u = DDBUFF(PSI,dvpdrdr)+Two_Over_R(r)*buffer(PSI,dvpdr) @@ -684,7 +720,9 @@ Subroutine Viscous_Force(buffer) ! r-direction; fluctuating If (compute_quantity(viscous_pforce_r) .or. & - compute_quantity(visc_work_pp)) Then + compute_quantity(visc_work_pp) .or. & + compute_quantity(curl_viscous_pforce_theta) .or. & + compute_quantity(curl_viscous_pforce_phi)) Then DO_PSI ! first, compute all the terms multiplied by mu @@ -715,7 +753,9 @@ Subroutine Viscous_Force(buffer) !Theta-direction; Fluctuating If (compute_quantity(viscous_pforce_theta) .or. & - compute_quantity(visc_work_pp)) Then + compute_quantity(visc_work_pp) .or. & + compute_quantity(curl_viscous_pforce_r) .or. & + compute_quantity(curl_viscous_pforce_phi)) Then DO_PSI ! first, compute all the terms multiplied by mu @@ -741,7 +781,9 @@ Subroutine Viscous_Force(buffer) !Phi-direction (fluctuating) If (compute_quantity(viscous_pforce_phi) .or. & - compute_quantity(visc_work_pp)) Then + compute_quantity(visc_work_pp) .or. & + compute_quantity(curl_viscous_pforce_r) .or. & + compute_quantity(curl_viscous_pforce_theta)) Then DO_PSI del2u = d2_fbuffer(PSI,dvpdrdr)+Two_Over_R(r)*fbuffer(PSI,dvpdr) @@ -770,7 +812,9 @@ Subroutine Viscous_Force(buffer) ! r-direction; mean If (compute_quantity(viscous_mforce_r) .or. & - compute_quantity(visc_work_mm)) Then + compute_quantity(visc_work_mm) .or. & + compute_quantity(curl_viscous_mforce_theta) .or. & + compute_quantity(curl_viscous_mforce_phi)) Then DO_PSI ! first, compute all the terms multiplied by mu @@ -799,7 +843,9 @@ Subroutine Viscous_Force(buffer) !Theta-direction; Mean If (compute_quantity(viscous_mforce_theta) .or. & - compute_quantity(visc_work_mm)) Then + compute_quantity(visc_work_mm) .or. & + compute_quantity(curl_viscous_mforce_r) .or. & + compute_quantity(curl_viscous_mforce_phi)) Then DO_PSI ! first, compute all the terms multiplied by mu @@ -826,7 +872,9 @@ Subroutine Viscous_Force(buffer) !Phi-direction (mean) If (compute_quantity(viscous_mforce_phi) .or. & compute_quantity(samom_diffusion) .or. & - compute_quantity(visc_work_mm)) Then + compute_quantity(visc_work_mm) .or. & + compute_quantity(curl_viscous_mforce_r) .or. & + compute_quantity(curl_viscous_mforce_theta)) Then DO_PSI del2u = d2_m0(PSI2,dvpdrdr)+Two_Over_R(r)*m0_values(PSI2,dvpdr) @@ -856,11 +904,4 @@ Subroutine Viscous_Force(buffer) End Subroutine Viscous_Force - Subroutine Finalize_Viscous_Force() - Implicit None - If (nvf .gt. 0) Then - DeAllocate(vforce_buffer) - Endif - End Subroutine Finalize_Viscous_Force - End Module Diagnostics_Linear_Forces 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/src/Diagnostics/curl_momentum_equation_codes.F b/src/Diagnostics/curl_momentum_equation_codes.F new file mode 100644 index 000000000..d8452215f --- /dev/null +++ b/src/Diagnostics/curl_momentum_equation_codes.F @@ -0,0 +1,146 @@ +! +! 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 :: 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(\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: $\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: $\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: $\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: $\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 + ! Note: the pressure gradient diagnostic codes are above + + 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 :: 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 :: 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 :: 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 :: 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_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+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+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+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+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 /////////////////////////////// + ! 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 :: 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+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+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+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+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+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$ 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 \ 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 +/