diff --git a/echemfem/solver.py b/echemfem/solver.py
index 29fee79..7dabd51 100644
--- a/echemfem/solver.py
+++ b/echemfem/solver.py
@@ -184,7 +184,9 @@ def __init__(
         self.poly_degreeU = self.poly_degree
         self.poly_degreep = self.poly_degree
         self.family = family
-        # DG penalization parameter. defined here for PMG
+        # DG penalization parameter
+        self.C_IP = Constant(10.)
+        # These DG parameters are defined here for PMG
         if p_penalty is None:
             p_penalty = Constant(p)
         self.penalty_degree = p_penalty  # Constant(self.poly_degree)
@@ -1487,8 +1489,7 @@ def mass_conservation_form(
         n = FacetNormal(mesh)
         if family == "DG":
             he = CellVolume(mesh) / FacetArea(mesh)
-            C_IP = 10.0
-            D_IP = C_IP * max_value(self.penalty_degree**2, 1) / he
+            D_IP = self.C_IP * max_value(self.penalty_degree**2, 1) / he
             K_IP = -1  # SIP
 
         a = 0.0
@@ -1762,8 +1763,7 @@ def gas_mass_conservation_form(self, X, test_fn, gas_params, u=None):
         n = FacetNormal(mesh)
         if family == "DG":
             he = CellVolume(mesh) / FacetArea(mesh)
-            C_IP = 10.0
-            D_IP = C_IP * max_value(self.penalty_degree**2, 1) / he
+            D_IP = self.C_IP * max_value(self.penalty_degree**2, 1) / he
             K_IP = -1  # SIP
 
         a = 0.0
@@ -1884,9 +1884,8 @@ def charge_conservation_form(self, u, test_fn, conc_params, W=None, i_bc=None):
         n = FacetNormal(mesh)
         if family == "DG":
             he = CellVolume(mesh) / FacetArea(mesh)
-            C_IP = Constant(10.0)
-            D_IP = C_IP * max_value(self.penalty_degree**2, 1) / he
-            D_IP2 = C_IP * max_value(self.penalty_degreeU**2, 1) / he
+            D_IP = self.C_IP * max_value(self.penalty_degree**2, 1) / he
+            D_IP2 = self.C_IP * max_value(self.penalty_degreeU**2, 1) / he
             K_IP = -1  # SIP
 
         a = 0.0
@@ -2067,9 +2066,8 @@ def potential_poisson_form(self, u, test_fn, conc_params, solid=False, W=None, i
         n = FacetNormal(mesh)
         if family == "DG":
             he = CellVolume(mesh) / FacetArea(mesh)
-            C_IP = 10.0
-            D_IP = C_IP * max_value(self.penalty_degree**2, 1) / he
-            D_IP2 = C_IP * max_value((self.penalty_degreeU)**2, 1) / he
+            D_IP = self.C_IP * max_value(self.penalty_degree**2, 1) / he
+            D_IP2 = self.C_IP * max_value((self.penalty_degreeU)**2, 1) / he
             K_IP = -1  # SIP
 
         if self.flow["porous"] and solid:
@@ -2255,10 +2253,9 @@ def liquid_pressure_form(self, p, test_fn, conc_params, u=None, W=None, i_bc=Non
         n = FacetNormal(mesh)
         if family == "DG":
             he = CellVolume(mesh) / FacetArea(mesh)
-            C_IP = 10.0
-            D_IP = C_IP * max_value(self.penalty_degree**2, 1) / he
-            D_IP2 = C_IP * max_value((self.penalty_degreep)**2, 1) / he
-            D_IP3 = C_IP * max_value((self.penalty_degreeU)**2, 1) / he
+            D_IP = self.C_IP * max_value(self.penalty_degree**2, 1) / he
+            D_IP2 = self.C_IP * max_value((self.penalty_degreep)**2, 1) / he
+            D_IP3 = self.C_IP * max_value((self.penalty_degreeU)**2, 1) / he
             K_IP = -1  # SIP
 
         a = 0.0
@@ -2405,8 +2402,7 @@ def gas_mass_conservation_form_pressure(
         n = FacetNormal(mesh)
         if family == "DG":
             he = CellVolume(mesh) / FacetArea(mesh)
-            C_IP = 10.0
-            D_IP2 = C_IP * max_value(self.penalty_degreep**2, 1) / he
+            D_IP2 = self.C_IP * max_value(self.penalty_degreep**2, 1) / he
             K_IP = -1  # SIP
 
         a = 0.0
@@ -2522,9 +2518,8 @@ def gas_pressure_form(self, p, test_fn, gas_params, u=None):  # currently unused
         n = FacetNormal(mesh)
         if family == "degree":
             he = CellVolume(mesh) / FacetArea(mesh)
-            C_IP = 10.0
-            D_IP = C_IP * max_value(self.penalty_degree**2, 1) / he
-            D_IP2 = C_IP * max_value((self.penalty_degreep)**2, 1) / he
+            D_IP = self.C_IP * max_value(self.penalty_degree**2, 1) / he
+            D_IP2 = self.C_IP * max_value((self.penalty_degreep)**2, 1) / he
             K_IP = -1  # SIP
 
         n_g = self.num_g