Merge pull request #113 from fastpm/ode

fix lightcone ode bug (seemingly from gsl 2.7)

libfastpm/cosmology.c

@@ -345,17 +345,28 @@ static ode_soln growth_ode_solve(double a, FastPMCosmology * c)
     yini[3] = 2 * yini[2];
 
     int status = gsl_odeiv2_driver_apply(drive, &aini, a, yini);
-    if (status != GSL_SUCCESS) {
-        fastpm_raise(-1, "Growth ODE unsuccesful at a=%g.", a);
-    }
-
     gsl_odeiv2_driver_free(drive);
 
     ode_soln soln;
-    soln.y0 = yini[0];
-    soln.y1 = yini[1];
-    soln.y2 = yini[2];
-    soln.y3 = yini[3];
+    if (status != GSL_SUCCESS) {
+        // if the ode solver failed and a >= aini, then there's a serious problem.
+        // if the ode solver failed and a < aini, this would either happen:
+        // 1) for the horizon interpolation table which is unused for a < aini anyway, so just return 0.
+        // 2) if the user input a < aini. In this case the results of the sim would be incorrect (FIXME).
+        if (a >= aini) {
+            fastpm_raise(-1, "Growth ODE unsuccesful at a=%g.", a);
+        } else {
+            soln.y0 = 0;
+            soln.y1 = 0;
+            soln.y2 = 0;
+            soln.y3 = 0;
+        }
+    } else {
+        soln.y0 = yini[0];
+        soln.y1 = yini[1];
+        soln.y2 = yini[2];
+        soln.y3 = yini[3];
+    }
 
     return soln;
 }

src/lua-runtime-fastpm.lua

@@ -78,32 +78,35 @@ function schema.omega_m.action (value)
 end
 
 -- check for bad input
-function schema.T_cmb.action (T_cmb)
-    if T_cmb ~= 0 then
+function schema.time_step.action (time_step)
+    function schema.T_cmb.action (T_cmb)
         function schema.growth_mode.action (growth_mode)
-            if growth_mode ~= 'ODE' then
+            if T_cmb ~= 0 and growth_mode ~= 'ODE' then
                 error("For a run with radiation (T_cmb > 0) use growth_mode='ODE' for accurate results.")
             end
+            -- enforce a<0.00625 (hard coded) when using ODE growth mode
+            if time_step[1] < 0.00625 and growth_mode == 'ODE' then
+                error("Cannot start the simulation at a<0.00625 when growth_mode=='ODE'.")
+            end
         end
-    end
-
-    function schema.m_ncdm.action (m_ncdm)
-        if #m_ncdm ~= 0 then
-            for i=2, #m_ncdm do
-                if m_ncdm[i] > m_ncdm[1] then
-                    error("Please input the heaviest ncdm particle first.")
+        function schema.m_ncdm.action (m_ncdm)
+            if #m_ncdm ~= 0 then
+                for i=2, #m_ncdm do
+                    if m_ncdm[i] > m_ncdm[1] then
+                        error("Please input the heaviest ncdm particle first.")
+                    end
                 end
-            end
-            function schema.n_shell.action (n_shell)
-                function schema.ncdm_freestreaming.action (ncdm_freestreaming)
-                    if ncdm_freestreaming and n_shell ~= 0 then
-                         error("For free-streaming ncdm use n_shell = 0 to turn off ncdm particles.")
+                function schema.n_shell.action (n_shell)
+                    function schema.ncdm_freestreaming.action (ncdm_freestreaming)
+                        if ncdm_freestreaming and n_shell ~= 0 then
+                             error("For free-streaming ncdm use n_shell = 0 to turn off ncdm particles.")
+                        end
                     end
                 end
-            end
-            function schema.ncdm_matterlike.action (ncdm_matterlike)
-                if not ncdm_matterlike and T_cmb == 0 then
-                     error("For a run with exact Omega_ncdm, T_cmb > 0 is required.")
+                function schema.ncdm_matterlike.action (ncdm_matterlike)
+                    if not ncdm_matterlike and T_cmb == 0 then
+                         error("For a run with exact Omega_ncdm, T_cmb > 0 is required.")
+                    end
                 end
             end
         end