From 226998e80efbd89b9319eac70bc42a44a17f0d33 Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Thu, 6 Jan 2022 16:46:36 -0300
Subject: [PATCH] SOURCES: Add source terms related to the MHD KEPES solver.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
 sources/sources.F90 | 59 ++++++++++++++++++++++++++++++---------------
 1 file changed, 40 insertions(+), 19 deletions(-)

diff --git a/sources/sources.F90 b/sources/sources.F90
index 8976909..9fa7399 100644
--- a/sources/sources.F90
+++ b/sources/sources.F90
@@ -121,6 +121,9 @@ module sources
     case("heglm", "HEGLM")
       glm_type = 2
       glm_name = "HEGLM"
+    case("kepes", "KEPES")
+      glm_type = 3
+      glm_name = "KEPES"
     case default
       glm_type = 0
       glm_name = "none"
@@ -556,7 +559,7 @@ module sources
       dh(2) = ady(pdata%meta%level)
       dh(3) = adz(pdata%meta%level)
 
-! add the EGLM-MHD source terms
+! add the EGLM-MHD and KEPES source terms
 !
       if (glm_type > 0) then
 
@@ -571,23 +574,7 @@ module sources
         du(imy,:,:,:) = du(imy,:,:,:) - db(:,:,:) * pdata%q(iby,:,:,:)
         du(imz,:,:,:) = du(imz,:,:,:) - db(:,:,:) * pdata%q(ibz,:,:,:)
 
-! update the energy equation
-!
-        if (ien > 0 .and. ibp > 0) then
-
-! calculate the gradient of divergence potential
-!
-          call gradient(dh(:), pdata%q(ibp,:,:,:), tmp(inx:inz,inx,:,:,:))
-
-! add the divergence potential source term to the energy equation, i.e.
-!     d/dt E + ∇.F = - B.(∇ψ)
-!
-          du(ien,:,:,:) = du(ien,:,:,:)                                        &
-                     - sum(pdata%q(ibx:ibz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)
-
-        end if ! ien > 0
-
-! add the HEGLM-MHD source terms
+! add the HEGLM-MHD and KEPES source terms
 !
         if (glm_type > 1) then
 
@@ -616,7 +603,41 @@ module sources
 
         end if ! glm_type > 1
 
-      end if ! glmtype > 0
+        if  (ibp > 0 .and. (ien > 0 .or. glm_type == 3)) then
+
+! calculate the gradient of divergence correcting field
+!
+          call gradient(dh(:), pdata%q(ibp,:,:,:), tmp(inx:inz,inx,:,:,:))
+
+          db(:,:,:) = sum(pdata%q(ivx:ivz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)
+
+! update the divergence correcting field
+!     d/dt ψ+ ∇.F = - (v.∇)ψ
+!
+          du(ibp,:,:,:) = du(ibp,:,:,:) - db(:,:,:)
+
+! update the energy equation
+!
+          if (ien > 0) then
+            if  (glm_type == 3) then
+
+! add the divergence potential source term to the energy equation, i.e.
+!     d/dt E + ∇.F = - ψ(v.∇)ψ
+!
+              du(ien,:,:,:) = du(ien,:,:,:) - pdata%q(ibp,:,:,:) * db(:,:,:)
+
+            else
+
+! add the divergence potential source term to the energy equation, i.e.
+!     d/dt E + ∇.F = - B.(∇ψ)
+!
+              du(ien,:,:,:) = du(ien,:,:,:) &
+                       - sum(pdata%q(ibx:ibz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)
+
+            end if
+          end if ! ien > 0
+        end if ! glm_type == 3
+      end if ! glm_type > 0
 
 ! if anomalous resistivity is enabled
 !