From 7cd2414faa117e427a28fd2a787aa48491c8b65e Mon Sep 17 00:00:00 2001
From: CGH0S7 <776459475@qq.com>
Date: Thu, 9 Apr 2026 19:17:39 +0800
Subject: [PATCH] Move constraint recomputation onto GPU path

---
 AMSS_NCKU_source/bssn_class.C | 60 +++++++++++++++++++++--------------
 1 file changed, 36 insertions(+), 24 deletions(-)

diff --git a/AMSS_NCKU_source/bssn_class.C b/AMSS_NCKU_source/bssn_class.C
index 9e085be..6733383 100644
--- a/AMSS_NCKU_source/bssn_class.C
+++ b/AMSS_NCKU_source/bssn_class.C
@@ -29,8 +29,12 @@ using namespace std;
 #include "rungekutta4_rout.h"
 #include "sommerfeld_rout.h"
 #include "getnp4.h"
-#include "shellfunctions.h"
-#include "parameters.h"
+#include "shellfunctions.h"
+#include "parameters.h"
+#ifdef USE_GPU
+#include "bssn_macro.h"
+#include "bssn_gpu.h"
+#endif
 
 #ifdef With_AHF
 #include "derivatives.h"
@@ -7282,12 +7286,15 @@ void bssn_class::Constraint_Out()
           MyList<Block> *BP = Pp->data->blb;
           while (BP)
           {
-            Block *cg = BP->data;
-            if (myrank == cg->rank)
-            {
-              f_compute_rhs_bssn(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2],
-                                 cg->fgfs[phi0->sgfn], cg->fgfs[trK0->sgfn],
-                                 cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], 
+            Block *cg = BP->data;
+            if (myrank == cg->rank)
+            {
+#ifdef USE_GPU
+              gpu_rhs(CALLED_BY_CONSTRAINT, myrank, RHS_PARA_CALLED_Constraint_Out);
+#else
+              f_compute_rhs_bssn(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2],
+                                 cg->fgfs[phi0->sgfn], cg->fgfs[trK0->sgfn],
+                                 cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], 
                                  cg->fgfs[gyy0->sgfn], cg->fgfs[gyz0->sgfn], cg->fgfs[gzz0->sgfn],
                                  cg->fgfs[Axx0->sgfn], cg->fgfs[Axy0->sgfn], cg->fgfs[Axz0->sgfn], 
                                  cg->fgfs[Ayy0->sgfn], cg->fgfs[Ayz0->sgfn], cg->fgfs[Azz0->sgfn],
@@ -7315,11 +7322,12 @@ void bssn_class::Constraint_Out()
                                  cg->fgfs[Gamzyy->sgfn], cg->fgfs[Gamzyz->sgfn], cg->fgfs[Gamzzz->sgfn],
                                  cg->fgfs[Rxx->sgfn], cg->fgfs[Rxy->sgfn], cg->fgfs[Rxz->sgfn], 
                                  cg->fgfs[Ryy->sgfn], cg->fgfs[Ryz->sgfn], cg->fgfs[Rzz->sgfn],
-                                 cg->fgfs[Cons_Ham->sgfn],
-                                 cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn],
-                                 cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn],
-                                 Symmetry, lev, ndeps, pre);
-            }
+                                 cg->fgfs[Cons_Ham->sgfn],
+                                 cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn],
+                                 cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn],
+                                 Symmetry, lev, ndeps, pre);
+#endif
+            }
             if (BP == Pp->data->ble)
               break;
             BP = BP->next;
@@ -7785,12 +7793,15 @@ void bssn_class::Interp_Constraint(bool infg)
           MyList<Block> *BP = Pp->data->blb;
           while (BP)
           {
-            Block *cg = BP->data;
-            if (myrank == cg->rank)
-            {
-              f_compute_rhs_bssn(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2],
-                                 cg->fgfs[phi0->sgfn], cg->fgfs[trK0->sgfn],
-                                 cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], 
+            Block *cg = BP->data;
+            if (myrank == cg->rank)
+            {
+#ifdef USE_GPU
+              gpu_rhs(CALLED_BY_CONSTRAINT, myrank, RHS_PARA_CALLED_Interp_Constraint);
+#else
+              f_compute_rhs_bssn(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2],
+                                 cg->fgfs[phi0->sgfn], cg->fgfs[trK0->sgfn],
+                                 cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], 
                                  cg->fgfs[gyy0->sgfn], cg->fgfs[gyz0->sgfn], cg->fgfs[gzz0->sgfn],
                                  cg->fgfs[Axx0->sgfn], cg->fgfs[Axy0->sgfn], cg->fgfs[Axz0->sgfn], 
                                  cg->fgfs[Ayy0->sgfn], cg->fgfs[Ayz0->sgfn], cg->fgfs[Azz0->sgfn],
@@ -7818,11 +7829,12 @@ void bssn_class::Interp_Constraint(bool infg)
                                  cg->fgfs[Gamzyy->sgfn], cg->fgfs[Gamzyz->sgfn], cg->fgfs[Gamzzz->sgfn],
                                  cg->fgfs[Rxx->sgfn], cg->fgfs[Rxy->sgfn], cg->fgfs[Rxz->sgfn], 
                                  cg->fgfs[Ryy->sgfn], cg->fgfs[Ryz->sgfn], cg->fgfs[Rzz->sgfn],
-                                 cg->fgfs[Cons_Ham->sgfn],
-                                 cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn],
-                                 cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn],
-                                 Symmetry, lev, ndeps, pre);
-            }
+                                 cg->fgfs[Cons_Ham->sgfn],
+                                 cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn],
+                                 cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn],
+                                 Symmetry, lev, ndeps, pre);
+#endif
+            }
             if (BP == Pp->data->ble)
               break;
             BP = BP->next;