ell_relax: make the time stepping factor configurable

And use different defaults for different FD stencils.

2 files changed, 16 insertions, 2 deletions

diff --git a/ell_relax.c b/ell_relax.c
index e56296e..313af4e 100644
--- a/ell_relax.c
+++ b/ell_relax.c
@@ -25,7 +25,10 @@
 #include "ell_relax.h"
 #include "log.h"
 
-#define CFL_FACT (1.0 / 7.0)
+static const double relax_factors[FD_STENCIL_MAX] = {
+    [0] = 1.0 / 5,
+    [1] = 1.0 / 7,
+};
 
 struct EllRelaxInternal {
     ptrdiff_t stride;
@@ -37,6 +40,7 @@ struct EllRelaxInternal {
     double *boundary_val_base[4];
 
     void (*derivatives_calc)(EllRelaxContext *ctx, double *dst, ptrdiff_t idx, ptrdiff_t stride);
+    double relax_factor;
 };
 
 static const struct {
@@ -237,7 +241,7 @@ static void boundaries_apply(EllRelaxContext *ctx)
 int mg2di_ell_relax_step(EllRelaxContext *ctx)
 {
     EllRelaxInternal *priv = ctx->priv;
-    const double cfl_fac = CFL_FACT * ctx->step[0] * ctx->step[1];
+    const double cfl_fac = priv->relax_factor * ctx->step[0] * ctx->step[1];
     int64_t start;
 
     start = gettime();
@@ -279,6 +283,11 @@ int mg2di_ell_relax_init(EllRelaxContext *ctx)
         return -EINVAL;
     }
 
+    if (ctx->relax_factor == 0.0)
+        priv->relax_factor = relax_factors[ctx->fd_stencil - 1];
+    else
+        priv->relax_factor = ctx->relax_factor;
+
     for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
         if (ctx->boundaries[i].type == ELL_RELAX_BC_TYPE_FIXDIFF) {
            for (int k = 0; k < ctx->domain_size[boundary_def[i].stride_idx]; k++)
diff --git a/ell_relax.h b/ell_relax.h
index 0e32015..4e97882 100644
--- a/ell_relax.h
+++ b/ell_relax.h
@@ -179,6 +179,11 @@ typedef struct EllRelaxContext {
     EllRelaxBoundary boundaries[4];
 
     /**
+     * The time stepping factor in relaxation.
+     */
+    double relax_factor;
+
+    /**
      * Values of the unknown function.
      *
      * Allocated by the solver in mg2di_ell_relax_alloc(), owned by the solver.