Use better initial guesses for solving.

1 files changed, 34 insertions, 1 deletions

diff --git a/src/maximal_slicing_axi_mg.c b/src/maximal_slicing_axi_mg.c
index 2601859..2020079 100644
--- a/src/maximal_slicing_axi_mg.c
+++ b/src/maximal_slicing_axi_mg.c
@@ -111,6 +111,7 @@ typedef struct MSMGContext {
     int64_t time_fine_boundaries;
     int64_t time_fine_mg2d;
     int64_t time_fine_export;
+    int64_t time_init_guess_eval;
 } MSMGContext;
 
 static MSMGContext *ms;
@@ -610,18 +611,40 @@ void msa_mg_eval(CCTK_ARGUMENTS)
         const int timestep = lrint(t * cctkGH->cctk_levfac[0] / ms->base_dt);
         int mol_substeps = MoLNumIntegratorSubsteps();
         int *mol_step = CCTK_VarDataPtr(cctkGH, 0, "MoL::MoL_Intermediate_Step");
+        int solver_idx;
 
         MG2DContext *solver;
 
         if (!mol_step || *mol_step <= 0)
             CCTK_WARN(0, "Invalid MoL step");
 
-        solver = cp->solver_eval[(cp->cur_step & 1) * mol_substeps + (mol_substeps - *mol_step)];
+        solver_idx = (cp->cur_step & 1) * mol_substeps + (mol_substeps - *mol_step);
+        solver = cp->solver_eval[solver_idx];
 
         start = gettime();
         fill_eq_coeffs(ms, cp, solver);
         ms->time_fine_fill_coeffs += gettime() - start;
 
+        {
+            start = gettime();
+
+            if (solver_idx) {
+                MG2DContext *solver_src = cp->solver_eval[solver_idx - 1];
+
+                for (int line = 0; line < solver->domain_size; line++) {
+                    memcpy(solver->u + line * solver->u_stride, solver_src->u + line * solver_src->u_stride,
+                           sizeof(*solver->u) * solver->domain_size);
+                }
+            } else {
+                for (int line = 0; line < solver->domain_size; line++)
+                    for (int i = 0; i < solver->domain_size; i++)
+                        solver->u[line * solver->u_stride + i] = lapse_mg[(cp->offset_left[1] + line) * cctk_lsh[0] + cp->offset_left[0] + i] - 1.0;
+            }
+
+        }
+
+        ms->time_init_guess_eval += gettime() - start;
+
         LOGDEBUG( "extrapolating BCs from t1=%g (%g) and t0=%g (%g)\n", lapse_prev1_time[reflevel], fact1, lapse_prev0_time[reflevel], fact0);
 
         start = gettime();
@@ -752,7 +775,17 @@ void msa_mg_solve(CCTK_ARGUMENTS)
                     lapse_mg[idx] = fact0 * lapse_prev0[idx] + fact1 * lapse_prev1[idx];
                 }
             }
+        } else {
+            /* use the solution from the coarser level as the initial guess */
+            CoordPatch *cp1 = get_coord_patch(ms, reflevel - 1);
+
+            ret = mg2d_init_guess(cp->solver, cp1->solver->u, cp1->solver->u_stride,
+                                  (size_t [2]){ cp1->solver->domain_size, cp1->solver->domain_size },
+                                  cp1->solver->step);
+            if (ret < 0)
+                CCTK_WARN(0, "Error setting the initial guess");
         }
+
         /* if the condition above was false, then lapse_mg should be filled by
          * prolongation from the coarser level
          * note that the reflection-boundary ghost points are not filled