Work in progress to optimize the iteration method

legendre_integration.c

@@ -33,6 +33,9 @@ double legendre_integration_A(gsl_spline2d* Spline_I2d, double tau, double u, do
 
     status = gsl_spline2d_eval_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_l);
 
+    //printf("u %lf tau %lf Ik_l %5.4e\n", u_prime[ii], tau, Ik_l);
+
+
     if (Ik_l < 0)
     {
       printf("Warning: interpolated Ik_l < 0 for tau=%5.4f u_prime=%5.4f\n", tau, u_prime[ii]);
@@ -43,7 +46,7 @@ double legendre_integration_A(gsl_spline2d* Spline_I2d, double tau, double u, do
     {
       printf("Error at line %d\n", __LINE__);
       printf("File %s\n", __FILE__);
-      printf("tau=%5.4f u=%5.4f\n", tau, u_prime[ii]);
+      printf("tau=%5.4f u=%5.4f  I^{k}(l)=%5.4e\n", tau, u_prime[ii], Ik_l);
 
       exit(1);
     }

main_program.c

@@ -213,7 +213,7 @@ int main(int argc, char* argv[])
       exit(1);
     }
 
-    solve_rte(kiter, seed);
+    optimize_taustep(kiter, seed);
   }
   else if (strcmp(method, "st85") == 0)
   {

solve_rte.c

 #include "iteration_functions.h"
+#include <gsl/gsl_interp2d.h>
 
 double I0_center_updown(double tau, double u);
 double I0_upward(double tau, double u);
@@ -11,12 +12,28 @@ double tau0;
 
 /*================================================================================*/
 
-void solve_rte(int k_iter, int seed_distribution)
+void optimize_taustep(int k_iter, int seed_distribution)
 {
+  Nstep_tau = 100;
+
+  tau0 = disktau / 2.;
+
+  while (solve_rte(k_iter, seed_distribution, Nstep_tau))
+  {
+    Nstep_tau = (int)Nstep_tau * 1.2;
+  }
+}
+
+/*================================================================================*/
+
+int solve_rte(int k_iter, int seed_distribution, int Nstep_tau)
+{
+
+  printf("\n\nReceveing Nstep_tau %d\n", Nstep_tau);
   int ii, jj, kk, status;
 
   Nstep_mu = 30;
-  Nstep_tau = 100;
+  /* Nstep_tau = 100;*/
 
   tau0 = disktau / 2.;
 
@@ -47,7 +64,6 @@ void solve_rte(int k_iter, int seed_distribution)
   {
     array_mu[ii - 1] = roots[ii];
     array_weights[ii - 1] = weights[ii];
-
     // printf("%lf 1 \n", array_mu[ii - 1]);
   }
 
@@ -109,6 +125,9 @@ void solve_rte(int k_iter, int seed_distribution)
       }
     }
   }
+
+
+
   /*=================================================================*/
   /*Setup for 2D interpolation*/
   /*=================================================================*/
@@ -129,10 +148,8 @@ void solve_rte(int k_iter, int seed_distribution)
   xacc_2d = gsl_interp_accel_alloc();
   yacc_2d = gsl_interp_accel_alloc();
 
-  step_tau = (2 * tau0) / (Nstep_tau - 1);
+  step_tau = (2 * tau0) / Nstep_tau;
 
-  // printf("Step tau %5.4f\n", step_tau);
-  /*double check_val=0;*/
 
   Iklr_intensity* ptr_data = malloc(sizeof(Iklr_intensity));
   Iklr_intensity* params;
@@ -173,13 +190,23 @@ void solve_rte(int k_iter, int seed_distribution)
     ptr_data->Spline_I2d_r_upstream = Spline_I2d_r_upstream;
     ptr_data->Spline_I2d_r_downstream = Spline_I2d_r_downstream;
 
+
+
+
+    printf("CHECK HERE !\n");
+    double zio=0;
+        int status = gsl_spline2d_eval_extrap_e(Spline_I2d_l_downstream, 2*tau0, 0.5, xacc_2d, yacc_2d, &zio);
+
+       printf("================> status %d val %lf\n", status, zio);
+
+
     /*=======================================================================*/
     /*Solve the RTE at the quadrature points*/
     /*=======================================================================*/
 
     status = GSL_SUCCESS + 1;
 
-    gsl_set_error_handler_off();
+    //gsl_set_error_handler_off();
 
     for (jj = 0; jj < Nstep_mu; jj++)
     {
@@ -204,53 +231,26 @@ void solve_rte(int k_iter, int seed_distribution)
       y[2] = 0;
       y[3] = 0;
 
-      /*   h_step=1e-7;
-         gsl_odeiv_evolve * e_rte       = gsl_odeiv_evolve_alloc (4);
-         gsl_odeiv_control * c_rte      = gsl_odeiv_control_y_new (1e-8, 1e-8);
-         gsl_odeiv_step * s_rte         = gsl_odeiv_step_alloc (T, 4);
-   */
       gsl_odeiv2_system sys = {RTE_Equations, jac, 4, params};
-      gsl_odeiv2_driver* d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk4, step_tau, 1e-4, 1e-4);
+      gsl_odeiv2_driver* d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk4, step_tau, 1e-7, 1e-7);
 
       while (tau < 2 * tau0)
       {
-        status = gsl_odeiv2_driver_apply_fixed_step(d, &tau, step_tau, 5, y);
-
-        while (status != GSL_SUCCESS)
-        {
-          tau = 0;
-          step_tau = step_tau * 0.75;
-          printf("step_tau now is %lf\n", step_tau);
-
-          // gsl_odeiv2_driver_free(d);
-
-          gsl_odeiv2_driver* d =
-              gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd, step_tau, 1e-3, 1e-3);
+        // printf("tau prima %lf\n", tau);
 
         status = gsl_odeiv2_driver_apply_fixed_step(d, &tau, step_tau, 5, y);
-        }
-
-        idx = set_array_idx(tau_prev, tau, array_tau, Nstep_tau);
 
         if (status)
         {
           printf("Error during integration: %s\n", gsl_strerror(status));
-          printf("Status %d\n", status);
-          exit(1);
+          // printf("Status %d\n", status);
+          return (1);
         }
 
-        // printf("tau_prev=%lf tau=%lf idx %d\n", tau_prev, tau, idx);
-
-        // status = gsl_odeiv_evolve_apply(e_rte, c_rte, s_rte, &sys, &tau, 2*tau0, &h_step, y);
-
-        // printf("tau=%lf  y[0]=%4.3e  y[1]=%4.3e y[2]=%4.3e y[3]=%4.3e idx %d\n", tau, y[0], y[1],
-        // y[2], y[3], idx);
+        idx = set_array_idx(tau_prev, tau, array_tau, Nstep_tau);
 
         if (idx > 0)
         {
-          // printf("tau=%lf  y[0]=%4.3e  y[1]=%4.3e y[2]=%4.3e y[3]=%4.3e idx %d\n", tau, y[0], y[1],
-          //      y[2], y[3], idx);
-
           ptr_Iklr[kk].Il_matrix_upstream[idx][jj] = y[0];
           ptr_Iklr[kk].Ir_matrix_upstream[idx][jj] = y[1];
 
@@ -263,10 +263,6 @@ void solve_rte(int k_iter, int seed_distribution)
 
       gsl_odeiv2_driver_free(d);
 
-      // gsl_odeiv_evolve_free (e_rte);
-      // gsl_odeiv_control_free (c_rte);
-      // gsl_odeiv_step_free (s_rte);
-
     } // End of loop over array_mu
 
   } // End of loop over k_iter
@@ -276,7 +272,20 @@ void solve_rte(int k_iter, int seed_distribution)
 
   compute_results(1, k_iter, Nstep_tau, Nstep_mu, array_mu, array_weights, ptr_Iklr, ptr_a, ptr_b);
 
-  exit(1);
+  free(Il_upstream_za);
+  free(Il_downstream_za);
+  free(Ir_upstream_za);
+  free(Ir_downstream_za);
+  free(array_tau);
+  free(ptr_Iklr);
+
+  gsl_spline2d_free(Spline_I2d_l_upstream);
+  gsl_spline2d_free(Spline_I2d_l_downstream);
+  gsl_spline2d_free(Spline_I2d_r_upstream);
+  gsl_spline2d_free(Spline_I2d_r_downstream);
+
+
+  return 0;
 }
 
 /*================================================================================*/
@@ -295,15 +304,14 @@ int RTE_Equations(double s, const double y[], double f[], void* params)
   gsl_spline2d* Spline_I2d_r_upstream = ((Iklr_intensity*)params)->Spline_I2d_r_upstream;
   gsl_spline2d* Spline_I2d_r_downstream = ((Iklr_intensity*)params)->Spline_I2d_r_downstream;
 
-  if (s > 2 * tau0 - 1e-5)
-  {
-    s = 2 * tau0;
-  }
 
   /*==============================================*/
   /*Equation for I_l  upstream*/
   /*==============================================*/
 
+  //printf("UNO  s=%10.7f\n", 2*tau0);
+
+
   f[0] = -1 / u * y[0] +
          3 / (8 * u) *
              (legendre_integration_A(Spline_I2d_l_upstream, s, u, array_u, weights_u) +
@@ -312,10 +320,16 @@ int RTE_Equations(double s, const double y[], double f[], void* params)
               legendre_integration_B(Spline_I2d_r_upstream, s, u, array_u, weights_u) +
               legendre_integration_B(Spline_I2d_r_downstream, 2 * tau0 - s, u, array_u, weights_u));
 
+
+
+
+
   /*==============================================*/
   /*Equation for I_r upstream*/
   /*==============================================*/
 
+
+
   f[1] = -1 / u * y[1] +
          3 / (8 * u) *
              (legendre_integration_C(Spline_I2d_l_upstream, s, u, array_u, weights_u) +
@@ -327,6 +341,8 @@ int RTE_Equations(double s, const double y[], double f[], void* params)
   /*Equation for I_l downstream*/
   /*==============================================*/
 
+
+
   f[2] = -1 / u * y[2] +
          3 / (8 * u) *
              (legendre_integration_A(Spline_I2d_l_downstream, s, u, array_u, weights_u) +
@@ -339,6 +355,8 @@ int RTE_Equations(double s, const double y[], double f[], void* params)
   /*Equation for I_r downstream*/
   /*==============================================*/
 
+
+
   f[3] = -1 / u * y[3] +
          3 / (8 * u) *
              (legendre_integration_C(Spline_I2d_l_downstream, s, u, array_u, weights_u) +