add CAHD

diff --git a/CodeGen/bssn.py b/CodeGen/bssn.py
index 082f3b8..1a62118 100644
--- a/CodeGen/bssn.py
+++ b/CodeGen/bssn.py
@@ -55,7 +55,14 @@ ad = dendro.ad
 kod = dendro.kod
 d2 = dendro.d2

-t = symbols("t")
+t = symbols("t") # time
+
+# added symbols used for CAHD
+ham = symbols('ham') # hamiltonian constraint violation
+C_CAHD = symbols('BSSN_CAHD_C') # coefficient for CAHD strength
+dt = symbols('dt') # simulation time step
+dx_i = symbols('dx_i') # spatial resolution of current grid
+dx_min = symbols('dx_min') # spatial resolution of finest grid

 dendro.set_metric(gt)
 igt = dendro.get_inverse_metric()
@@ -68,7 +75,7 @@ eta_func = (
 )


-def bssn_puncture_gauge(eta_damp, isStaged=False, prefix="", sslGaugeCondition=False):
+def bssn_puncture_gauge(eta_damp, isStaged=False, prefix="", sslGaugeCondition=False, enableCAHD=False):
     """
     BSSN puncture gauge (HAD/ traditional BSSN puncture gaugue) with const eta damping
     """
@@ -100,6 +107,10 @@ def bssn_puncture_gauge(eta_damp, isStaged=False, prefix="", sslGaugeCondition=F
         gt_rhs = dendro.lie(b, gt, weight) - 2 * a * At

         chi_rhs = dendro.lie(b, chi, weight) + Rational(2, 3) * (chi * a * K)
+
+        if enableCAHD:
+          # turn on curvature-adjusted Hamiltonian-constraint damping
+          chi_rhs += C_CAHD * chi**2 * (dt * dx_i/dx_min) * ham

         AikAkj = Matrix(
             [
@@ -261,6 +272,10 @@ def bssn_puncture_gauge(eta_damp, isStaged=False, prefix="", sslGaugeCondition=F
         gt_rhs = dendro.lie(b, gt, weight) - 2 * a * At

         chi_rhs = dendro.lie(b, chi, weight) + Rational(2, 3) * (chi * a * K)
+
+        if enableCAHD:
+          # turn on curvature-adjusted Hamiltonian-constraint damping
+          chi_rhs += C_CAHD * chi**2 * (dt * dx_i/dx_min) * ham

         AikAkj = Matrix(
             [
@@ -701,11 +716,14 @@ def bssn_rochester_puncture_gauge(
             dendro.generate_separate([outs[i]], [vnames[i]], "[pp]")


-def main(staged_type, gauge, eta_damp, prefix, enable_ssl):
+def main(staged_type, gauge, eta_damp, prefix, enable_ssl, enable_cahd):

     if enable_ssl:
         print("// CODEGEN: SSL was enabled, adding term to gauge condition!")

+    if enable_cahd:
+        print("// CODEGEN: CAHD was enabled, adding damping term to chi!")
+
     if staged_type == "staged":
         print("//Codgen: generating staged version ")
         if gauge == "rochester":
@@ -721,10 +739,10 @@ def main(staged_type, gauge, eta_damp, prefix, enable_ssl):
             print("//Codgen: using standard gauge")
             if eta_damp == "func":
                 print("//Codgen: using eta func damping")
-                bssn_puncture_gauge(eta_func, True, prefix, enable_ssl)
+                bssn_puncture_gauge(eta_func, True, prefix, enable_ssl, enable_cahd)
             else:
                 print("//Codgen: using eta const damping")
-                bssn_puncture_gauge(eta, True, prefix, enable_ssl)
+                bssn_puncture_gauge(eta, True, prefix, enable_ssl, enable_cahd)

     else:
         print("//Codgen: generating unstage version ")
@@ -741,10 +759,10 @@ def main(staged_type, gauge, eta_damp, prefix, enable_ssl):
             print("//Codgen: using standard gauge")
             if eta_damp == "func":
                 print("//Codgen: using eta func damping")
-                bssn_puncture_gauge(eta_func, False, prefix, enable_ssl)
+                bssn_puncture_gauge(eta_func, False, prefix, enable_ssl, enable_cahd)
             else:
                 print("//Codgen: using eta const damping")
-                bssn_puncture_gauge(eta, False, prefix, enable_ssl)
+                bssn_puncture_gauge(eta, False, prefix, enable_ssl, enable_cahd)


 if __name__ == "__main__":
@@ -787,6 +805,12 @@ if __name__ == "__main__":
         action="store_true",
         help="Whether or not to generate with the SSL code",
     )
+    parser.add_argument(
+        "-c",
+        "--enable_cahd",
+        action="store_true",
+        help="Whether or not to enable CAHD",
+    )

     args = parser.parse_args()