From 3def71f6ed726c8543c4b5d85ad096036f9f509d Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Mon, 4 Apr 2022 08:52:10 -0300
Subject: [PATCH] USER_PROBLEMS: Add reconnection rate terms through YZ
 boundaries.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
 sources/user_problem.F90 | 100 +++++++++++++++++++++++++++++++++++++--
 1 file changed, 97 insertions(+), 3 deletions(-)

diff --git a/sources/user_problem.F90 b/sources/user_problem.F90
index 2bbc8c9..bb8e95d 100644
--- a/sources/user_problem.F90
+++ b/sources/user_problem.F90
@@ -1531,7 +1531,8 @@ module user_problem
   subroutine user_statistics(time)
 
     use blocks     , only : block_meta, block_data, list_data
-    use coordinates, only : nn => bcells, nb, ne
+    use coordinates, only : nn => bcells, nb, ne, nbm, nep, periodic
+    use coordinates, only : xmin, xmax
     use coordinates, only : adx, ady, adz, advol, ay, yarea
     use equations  , only : ivx, ivy, ivz, ibx, iby, ibz, ibp
     use helpers    , only : print_message, flush_and_sync
@@ -1554,7 +1555,7 @@ module user_problem
 !$  integer :: omp_get_thread_num
 
     integer                     :: ni, nl, nu, nm, np, status
-    real(kind=8)                :: dvol, dxz
+    real(kind=8)                :: dvol, dyz, dxz
     real(kind=8), dimension(3)  :: dh
     real(kind=8), dimension(nn) :: yc
     real(kind=8), dimension(32) :: rterms
@@ -1604,6 +1605,7 @@ module user_problem
 
       dvol  = advol(pmeta%level)
 
+      dyz   = ady(pmeta%level) * adz(pmeta%level)
       dxz   = adx(pmeta%level) * adz(pmeta%level)
 
       yc(:) = pmeta%ymin + ay(pmeta%level,:)
@@ -1641,6 +1643,98 @@ module user_problem
                    + sum(pdata%u(ibx:ibz,nb:ne,nl:nu,  :  )**2) * dvol
 #endif /* NDIMS == 3 */
 
+        if (.not. periodic(1)) then
+          if (pmeta%xmin <= xmin) then
+
+! advection of magnetic energy across the lower X-boundary
+!
+            rterms(14) = rterms(14) &
+#if NDIMS == 3
+                     + sum(sum(pdata%q(ibx:ibz,nbm:nb,nl:nu,nb:ne)**2,1) &
+                             * pdata%q(ivx    ,nbm:nb,nl:nu,nb:ne)) * dyz
+#else /* NDIMS == 3 */
+                     + sum(sum(pdata%q(ibx:ibz,nbm:nb,nl:nu,  :  )**2,1) &
+                             * pdata%q(ivx    ,nbm:nb,nl:nu,  :  )) * dyz
+#endif /* NDIMS == 3 */
+
+            rterms(17) = rterms(17) &
+#if NDIMS == 3
+                     - sum(sum(pdata%q(ivx:ivz,nbm:nb,nl:nu,nb:ne) &
+                             * pdata%q(ibx:ibz,nbm:nb,nl:nu,nb:ne),1) &
+                             * pdata%q(ibx    ,nbm:nb,nl:nu,nb:ne)) * dyz
+#else /* NDIMS == 3 */
+                     - sum(sum(pdata%q(ivx:ivz,nbm:nb,nl:nu,  :  ) &
+                             * pdata%q(ibx:ibz,nbm:nb,nl:nu,  :  ),1) &
+                             * pdata%q(ibx    ,nbm:nb,nl:nu,  :  )) * dyz
+#endif /* NDIMS == 3 */
+
+          if (resistive) then
+
+! diffusion of magnetic energy across the lower X-boundary
+!
+            rterms(20) = rterms(20) &
+#if NDIMS == 3
+                         - sum(cur(2,nbm:nb,nl:nu,nb:ne) &
+                       * pdata%q(ibz,nbm:nb,nl:nu,nb:ne) &
+                             - cur(3,nbm:nb,nl:nu,nb:ne) &
+                       * pdata%q(iby,nbm:nb,nl:nu,nb:ne)) * dyz
+#else /* NDIMS == 3 */
+                         - sum(cur(2,nbm:nb,nl:nu,  :  ) &
+                       * pdata%q(ibz,nbm:nb,nl:nu,  :  ) &
+                             - cur(3,nbm:nb,nl:nu,  :  ) &
+                       * pdata%q(iby,nbm:nb,nl:nu,  :  )) * dyz
+#endif /* NDIMS == 3 */
+
+            end if ! resistivity
+
+          end if
+
+          if (pmeta%xmax >= xmax) then
+
+! advection of magnetic energy across the upper X-boundary
+!
+            rterms(14) = rterms(14) &
+#if NDIMS == 3
+                     - sum(sum(pdata%q(ibx:ibz,ne:nep,nl:nu,nb:ne)**2,1) &
+                             * pdata%q(ivx    ,ne:nep,nl:nu,nb:ne)) * dyz
+#else /* NDIMS == 3 */
+                     - sum(sum(pdata%q(ibx:ibz,ne:nep,nl:nu,  :  )**2,1) &
+                             * pdata%q(ivx    ,ne:nep,nl:nu,  :  )) * dyz
+#endif /* NDIMS == 3 */
+
+            rterms(17) = rterms(17) &
+#if NDIMS == 3
+                     + sum(sum(pdata%q(ivx:ivz,ne:nep,nl:nu,nb:ne) &
+                             * pdata%q(ibx:ibz,ne:nep,nl:nu,nb:ne),1) &
+                             * pdata%q(ibx    ,ne:nep,nl:nu,nb:ne)) * dyz
+#else /* NDIMS == 3 */
+                     + sum(sum(pdata%q(ivx:ivz,ne:nep,nl:nu,  :  ) &
+                             * pdata%q(ibx:ibz,ne:nep,nl:nu,  :  ),1) &
+                             * pdata%q(ibx    ,ne:nep,nl:nu,  :  )) * dyz
+#endif /* NDIMS == 3 */
+
+          if (resistive) then
+
+! diffusion of magnetic energy across the upper X-boundary
+!
+            rterms(20) = rterms(20) &
+#if NDIMS == 3
+                         + sum(cur(2,ne:nep,nl:nu,nb:ne) &
+                       * pdata%q(ibz,ne:nep,nl:nu,nb:ne) &
+                             - cur(3,ne:nep,nl:nu,nb:ne) &
+                       * pdata%q(iby,ne:nep,nl:nu,nb:ne)) * dyz
+#else /* NDIMS == 3 */
+                         + sum(cur(2,ne:nep,nl:nu,  :  ) &
+                       * pdata%q(ibz,ne:nep,nl:nu,  :  ) &
+                             - cur(3,ne:nep,nl:nu,  :  ) &
+                       * pdata%q(iby,ne:nep,nl:nu,  :  )) * dyz
+#endif /* NDIMS == 3 */
+
+            end if ! resistivity
+
+          end if
+        end if
+
         if (pmeta%ymin <= ylo .and. ylo < pmeta%ymax) then
           ni = nl
           nm = ni - 1
@@ -1696,7 +1790,7 @@ module user_problem
                      + sum(pdata%q(ibx:ibz,nb:ne,nm:ni,  :  )**2) * dxz
 #endif /* NDIMS == 3 */
 
-! advection of magnetic energy
+! advection of magnetic energy across the lower Y-boundary
 !
           rterms(15) = rterms(15) &
 #if NDIMS == 3