Merge branch 'master' into profile_trapping to update to unified iterations

 * \brief Implementation of the calculation for a cluster of spheres.

 */

#include <chrono>

#include <cmath>

#include <cstdio>

#include <exception>

#include <fstream>

    int nsph = gconf->number_of_spheres;

    // Sanity check on number of sphere consistency, should always be verified

    if (s_nsph == nsph) {

      char virtual_line[256];

      sprintf(virtual_line, "%.5g.\n", sconf->get_particle_radius(gconf));

      message = "INFO: particle radius is " + (string)virtual_line;

      logger->log(message);

      ScatteringAngles *p_scattering_angles = new ScatteringAngles(gconf);

      double wp = sconf->wp;

      // ClusterOutputInfo : Thread 0 of MPI process 0 allocates the memory to

      logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n");

      time_logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n");

      str(sconf, cid->c1);

      str(sconf->_rcf, cid->c1);

      thdps(cid->c1->lm, cid->zpv);

      double exdc = sconf->exdc;

      double exri = sqrt(exdc);

      }

      // do the first iteration on jxi488 separately, since it seems to be different from the others

      int jxi488 = 1;

      // not anymore, now this iteration is at the same level as others

      int jxi488;

      int initialmaxrefiters = cid->maxrefiters;

      chrono::time_point<chrono::high_resolution_clock> start_iter_1 = chrono::high_resolution_clock::now();

#ifdef USE_NVTX

      nvtxRangePush("First iteration");

#endif

      //chrono::time_point<chrono::high_resolution_clock> start_iter_1 = chrono::high_resolution_clock::now();

// #ifdef USE_NVTX

//       nvtxRangePush("First iteration");

// #endif

      // use these pragmas, which should have no effect on parallelism, just to push OMP nested levels at the same level also in the first wavelength iteration

      int jer = 0;

#pragma omp parallel

      {

#pragma omp single

	{

	  jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, cid, p_output, output_path, vtppoanp);

	} // OMP sinlge

      } // OMP parallel

#ifdef USE_NVTX

      nvtxRangePop();

#endif

      chrono::time_point<chrono::high_resolution_clock> end_iter_1 = chrono::high_resolution_clock::now();

      elapsed = start_iter_1 - t_start;

      string message = "INFO: Calculation setup took " + to_string(elapsed.count()) + "s.\n";

      logger->log(message);

      time_logger->log(message);

      elapsed = end_iter_1 - start_iter_1;

      message = "INFO: First iteration took " + to_string(elapsed.count()) + "s.\n";

      logger->log(message);

      time_logger->log(message);

// #pragma omp parallel

//       {

// #pragma omp single

// 	{

// 	  jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, cid, p_output, output_path, vtppoanp);

// 	} // OMP single

//       } // OMP parallel

// #ifdef USE_NVTX

//       nvtxRangePop();

// #endif

      //chrono::time_point<chrono::high_resolution_clock> end_iter_1 = chrono::high_resolution_clock::now();

      //elapsed = start_iter_1 - t_start;

      // string message = "INFO: Calculation setup took " + to_string(elapsed.count()) + "s.\n";

      // logger->log(message);

      // time_logger->log(message);

      // elapsed = end_iter_1 - start_iter_1;

      // message = "INFO: First iteration took " + to_string(elapsed.count()) + "s.\n";

      // logger->log(message);

      // time_logger->log(message);

      /* for the next iterations, just always do maxiter iterations, assuming the accuracy is good enough */

      cid->refinemode = 0;

      // cid->refinemode = 0;

      /* add an extra iteration for margin, if this does not exceed initialmaxrefiters */

      // if (cid->maxrefiters < initialmaxrefiters) cid->maxrefiters++;

      if (jer != 0) {

	// First loop failed. Halt the calculation.

	fclose(timing_file);

	delete time_logger;

	delete p_output;

	delete p_scattering_angles;

	delete cid;

	delete logger;

	delete sconf;

	delete gconf;

	return;

      }

      // if (jer != 0) {

      // 	// First loop failed. Halt the calculation.

      // 	fclose(timing_file);

      // 	delete time_logger;

      // 	delete p_output;

      // 	delete p_scattering_angles;

      // 	delete cid;

      // 	delete logger;

      // 	delete sconf;

      // 	delete gconf;

      // 	return;

      // }

      //==================================================

      // do the first outputs here, so that I open here the new files, afterwards I only append

      //==================================================

      // How should we handle this, when first iteration is not treated specially anymore? This should be ok, just write what was put in vtppoanp on initialisation, even if no actual calc was done yet. This creates the file nonetheless, 

      vtppoanp->write_to_disk(output_path + "/c_TPPOAN");

      delete vtppoanp;

	}

#pragma omp barrier

	// ok, now I can actually start the parallel calculations

	for (int ixi488=2; ixi488<=cid_2->number_of_scales; ixi488 +=myMPIstride) {

	// we now start the loop from 1, so that the first iteration is run exactly as the others

	for (int ixi488=1; ixi488<=cid_2->number_of_scales; ixi488 +=myMPIstride) {

	  // the parallel loop over MPI processes covers a different set of indices for each thread

#pragma omp barrier

	  int myjxi488 = ixi488+myompthread;

	    // thread 0 writes its virtual files, now including contributions from all threads, to disk, and deletes them

	    // p_outarray[0]->append_to_disk(output_path + "/c_OCLU");

	    // delete p_outarray[0];

	    // ******************************************************

	    // if this is the very first time, we should actually use

	    // ->write_to_disk, not ->append_to_disk

	    // ******************************************************

	    vtppoanarray[0]->append_to_disk(output_path + "/c_TPPOAN");

	    delete vtppoanarray[0];

      // make sure all threads align here: I don't want the following loop to accidentally start for thread 0, possibly modifying some variables before they are copied by all other threads

#pragma omp barrier

      // ok, now I can actually start the parallel calculations

      for (int ixi488=2; ixi488<=cid_2->number_of_scales; ixi488 +=myMPIstride) {

      // we now start the loop from 1, so that the first iteration is treated exactly on equal footing as others

      for (int ixi488=1; ixi488<=cid_2->number_of_scales; ixi488 +=myMPIstride) {

	// the parallel loop over MPI processes covers a different set of indices for each thread

#pragma omp barrier

	int myjxi488 = ixi488 + myjxi488startoffset + myompthread;

#endif

}

int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, ClusterIterationData *cid, ClusterOutputInfo *output, const string& output_path, VirtualBinaryFile *vtppoanp)

{

int cluster_jxi488_cycle(

  int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf,

  ScatteringAngles *sa, ClusterIterationData *cid, ClusterOutputInfo *output,

  const string& output_path, VirtualBinaryFile *vtppoanp

) {

  int nxi = sconf->number_of_scales;

  const dcomplex cc0 = 0.0 + I * 0.0;

  const double pi = acos(-1.0);

  char virtual_line[256];

  string message = "INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n";

  Logger *logger = new Logger(LOG_DEBG);

  int jer = 0;

  int lcalc = 0;

  int jaw = 1;

  int li = gconf->li;

  int le = gconf->le;

  int lm = 0;

  if (le > lm) lm = le;

  if (li > lm) lm = li;

  int nsph = sconf->number_of_spheres;

  np_int mxndm = gconf->mxndm;

  int iavm = gconf->iavm;

  int inpol = gconf->in_pol;

  int npnt = gconf->npnt;

  int npntts = gconf->npntts;

  int isam = gconf->isam;

  int jwtm = gconf->jwtm;

  np_int ndit = 2 * nsph * cid->c1->nlim;

  // this is now part of cid, so don't mess with it here, just copy it by reference

  bool &is_first_scale = cid->is_first_scale;

  const int nsph = sconf->number_of_spheres;

  const np_int mxndm = gconf->mxndm;

  const int iavm = gconf->iavm;

  const int inpol = gconf->in_pol;

  const int npnt = gconf->npnt;

  const int npntts = gconf->npntts;

  const int isam = gconf->isam;

  const int jwtm = gconf->jwtm;

  int isq, ibf;

  int last_configuration;

  int num_configs = sconf->configurations;

  int ndirs = sa->nkks;

  const int num_configs = sconf->configurations;

  const int ndirs = sa->nkks;

  int oindex = -1;

  int jindex = jxi488 - output->first_xi + 1;

  double xi = sconf->get_scale(jxi488 - 1);

  double exdc = sconf->exdc;

  double exri = sqrt(exdc);

  int idfc = (int)sconf->idfc;

  const int idfc = (int)sconf->idfc;

  double vkarg = 0.0;

  if (idfc >= 0) {

    cid->vk = xi * cid->wn;

    output->vec_vk[jindex - 1] = cid->vk;

    output->vec_xi[jindex - 1] = xi;

  }

  // Dynamic order check

  const int max_li = gconf->li;

  const int max_le = gconf->le;

  const double alamb = 2.0 * pi / cid->vk;

  double size_par_li = 2.0 * pi * sqrt(exdc) * sconf->get_max_radius() / alamb;

  int recommended_li = 4 + (int)ceil(size_par_li + 4.05 * pow(size_par_li, 1.0 / 3.0));

  double size_par_le = 2.0 * pi * sqrt(exdc) * sconf->get_particle_radius(gconf) / alamb;

  int recommended_le = 1 + (int)ceil(size_par_le + 11.0 * pow(size_par_le, 1.0 / 3.0));

  if (recommended_li != cid->c1->li || recommended_le != cid->c1->le) {

    if (recommended_li > cid->c1->li) {

      message = "WARNING: internal order " + to_string(cid->c1->li) + " for scale iteration "

	+ to_string(jxi488) + " too low (recommended order is " + to_string(recommended_li)

	+ ").\n";

      logger->log(message, LOG_WARN);

    } else if (recommended_li < cid->c1->li) {

      message = "INFO: lowering internal order from " + to_string(cid->c1->li) + " to "

	+ to_string(recommended_li) + " for scale iteration " + to_string(jxi488) + ".\n";

      logger->log(message, LOG_INFO);

    }

    if (recommended_le > cid->c1->le) {

      message = "WARNING: external order " + to_string(cid->c1->le) + " for scale iteration "

	+ to_string(jxi488) + " too low (recommended order is " + to_string(recommended_le)

	+ ").\n";

      logger->log(message, LOG_WARN);

    } else if (recommended_le < cid->c1->le) {

      message = "INFO: lowering external order from " + to_string(cid->c1->le) + " to "

	+ to_string(recommended_le) + " for scale iteration " + to_string(jxi488) + ".\n";

      logger->log(message, LOG_INFO);

    }

    if (recommended_li < max_li || recommended_le < max_le) {

      int new_li = (recommended_li < max_li) ? recommended_li : max_li;

      int new_le = (recommended_le < max_le) ? recommended_le : max_le;

      cid->update_orders(sconf->_rcf, new_li, new_le);

      is_first_scale = true;

      jaw = 1;

      cid->refinemode = 2;

    }

  }

  int li = cid->c1->li;

  int le = cid->c1->le;

  int lm = cid->c1->lm;

  np_int ndit = 2 * nsph * cid->c1->nlim;

  // End of dynamic order check

  hjv(exri, vkarg, jer, lcalc, cid->arg, cid->c1);

  if (jer != 0) {

    output->vec_ier[jindex - 1] = 1;

	for (int ic = 0; ic < ici; ic++)

	  cid->c1->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, jxi488 - 1);

      } else {

	if (jxi488 == 1) {

	if (is_first_scale) {

	  for (int ic = 0; ic < ici; ic++)

	    cid->c1->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, 0);

	}

    }

  }

#endif

  cid->refinemode = 0;

#ifdef DEBUG_AM

  VirtualAsciiFile *outam2 = new VirtualAsciiFile();

  string outam2_name = output_path + "/c_AM2_JXI" + to_string(jxi488) + ".txt";

    double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0;

    for (int jph484 = 1; jph484 <= sa->nph; jph484++) {

      int jw = 0;

      if (sa->nk != 1 || jxi488 <= 1) {

      if (sa->nk != 1 || is_first_scale) {

	upvmp(th, ph, 0, cost, sint, cosp, sinp, cid->u, cid->upmp, cid->unmp);

	if (isam >= 0) {

	  wmamp(

	  raba(cid->c1->le, cid->c1->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);

	  jw = 1;

	}

      } else { // label 180, NK == 1 AND JXI488 == 1

      } else { // label 180, NK == 1 AND JXI488 > 1

	if (isam >= 0) {

	  // label 182

	  apc(cid->zpv, cid->c1->le, cid->c1->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);

	    if (phs > 360.0) phs -= 360.0;

	  }

	  // label 188

	  bool goto190 = (sa->nks == 1 && (jxi488 > 1 || jth486 > 1 || jph484 > 1));

	  bool goto190 = (sa->nks == 1 && (!(is_first_scale) || jth486 > 1 || jph484 > 1));

	  if (!goto190) {

	    upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, cid->us, cid->upsmp, cid->unsmp);

	    if (isam >= 0)

		    );

	  }

	  // label 190

	  if (sa->nkks != 1 || jxi488 <= 1) {

	  if (sa->nkks != 1 || is_first_scale) {

	    upvsp(

		  cid->u, cid->upmp, cid->unmp, cid->us, cid->upsmp, cid->unsmp, cid->up, cid->un, cid->ups, cid->uns,

		  cid->duk, isq, ibf, cid->scan, cid->cfmp, cid->sfmp, cid->cfsp, cid->sfsp

    } // jph484 loop

    th += sa->thstp;

  } // jth486 loop

  // at the end of this iteration make sure to set is_first_scale to false, the next iteration will reinitialise _only_ if the order changes

  is_first_scale = 0;

#ifdef USE_NVTX

  nvtxRangePop();

#endif

  proc_device = 0;

#endif

  // the first time the object is used, it will be a "first iteration", to ensure data structures are properly initialised.

  is_first_scale = 1;

  // In the first iteration, if refinement is enabled, determine the number of refinement iterations required to arrive at the target accuracy (if achievable in a reasonable number of iterations)

  refinemode = 2;

  // maxrefiters and accuracygoal should be configurable and preferably set somewhere else

  number_of_scales = rhs.number_of_scales;

  proc_device = rhs.proc_device;

  is_first_scale = rhs.is_first_scale;

  refinemode = rhs.refinemode;

  maxrefiters = rhs.maxrefiters;

  accuracygoal = rhs.accuracygoal;

  proc_device = 0;

#endif

  MPI_Bcast(&refinemode, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&is_first_scale, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);

  MPI_Bcast(&maxrefiters, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&accuracygoal, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

}

  MPI_Bcast(&xiblock, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&number_of_scales, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&refinemode, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&is_first_scale, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);

  MPI_Bcast(&maxrefiters, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&accuracygoal, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

}

  delete[] cmullr;

  delete[] cmul;

}

int ClusterIterationData::update_orders(double **rcf, int inner_order, int outer_order) {

  int result = 0;

  int old_lm = c1->lm;

  np_int old_ndit = 2 * c1->nsph * c1->nlim;

  ((ParticleDescriptorCluster *)c1)->update_orders(inner_order, outer_order);

  const int ndi = c1->ndi;

  const np_int ndit = (np_int)c1->ndit;

  for (int zi = 0; zi < old_lm; zi++) {

    for (int zj = 0; zj < 3; zj++) {

      for (int zk = 0; zk < 2; zk++) {

	delete[] zpv[zi][zj][zk];

      }

      delete[] zpv[zi][zj];

    }

    delete[] zpv[zi];

  }

  delete[] zpv;

  zpv = new double***[c1->lm];

  for (int zi = 0; zi < c1->lm; zi++) {

    zpv[zi] = new double**[3];

    for (int zj = 0; zj < 3; zj++) {

      zpv[zi][zj] = new double*[2];

      for (int zk = 0; zk < 2; zk++) {

	zpv[zi][zj][zk] = new double[2]();

      }

    }

  }

  str(rcf, c1);

  thdps(c1->lm, zpv);

  delete[] am;

  delete[] am_vector;

  am_vector = new dcomplex[ndit * ndit]();

  am = new dcomplex*[ndit];

  for (int ai = 0; ai < ndit; ai++) {

    am[ai] = am_vector + ai * ndit;

  }

  return result;

}

// >>> END OF ClusterIterationData CLASS IMPLEMENTATION <<<

   * \return descriptor_type: `string` The descriptor type name.

   */

  std::string get_descriptor_type() override { return "cluster descriptor"; }

  /*! \brief Update the field expansion orders.

   *

   * \param inner_order: `int` The new inner expansion order to be set.

   * \param outer_order: `int` The new outer expansion order to be set.

   * \return result: `int` An exit code (0 if successful).

   */

  int update_orders(int inner_order, int outer_order);

};

/*! \brief The data structure describing a particle model for a sphere with inclusions.

   * \return descriptor_type: `string` The descriptor type name.

   */

  std::string get_descriptor_type() override { return "inclusion descriptor"; }

  /*! \brief Update the field expansion orders.

   *

   * \param inner_order: `int` The new inner expansion order to be set.

   * \param outer_order: `int` The new outer expansion order to be set.

   * \return result: `int` An exit code (0 if successful).

   */

  int update_orders(int inner_order, int outer_order);

};

/*! \brief The data structure describing a spherical particle model.

   * \return descriptor_type: `string` The descriptor type name.

   */

  std::string get_descriptor_type() override { return "sphere descriptor"; }

  /*! \brief Update the field expansion order.

   *

   * \param order: `int` The new field expansion order to be set.

   * \return result: `int` An exit code (0 if successful).

   */

  int update_order(int order);

};

/*! \brief A data structure representing the angles to be evaluated in the problem.

  dcomplex ***_dc0_matrix;

  //! \brief Vector of sphere radii expressed in m, with size [N_SPHERES].

  double *_radii_of_spheres;

  //! \brief Matrix of fractional transition radii with size [N_SPHERES x LAYERS].

  double **_rcf;

  //! \brief Vector of sphere ID numbers, with size [N_SPHERES].

  int *_iog_vec;

  //! \brief Vector of layer numbers for every sphere, with size [N_SPHERES].

  //! \brief Vector of layer numbers for every sphere, with size [CONFIGURATIONS].

  int *_nshl_vec;

  //! \brief Vector of scale parameters, with size [N_SCALES].

  double *_scale_vec;

  const int& max_layers = _max_layers;

  //! \brief Read only view on flag to control whether to add an external layer.

  const bool& use_external_sphere = _use_external_sphere;

  //! \brief Matrix of fractional transition radii with size [CONFIGURATIONS x LAYERS].

  double **_rcf;

  /*! \brief Build a scatterer configuration structure.

   *

   */

  int get_iog(int index) { return _iog_vec[index]; }

  /*! \brief Get the maximum radius of the sphere components.

   *

   * \return radius: `double` The radius of the largest sphere.

   */

  double get_max_radius();

  /*! \brief Get the number of layers for a given configuration.

   *

   * This is a specialized function to get the number of layers in a specific

   */

  int get_nshl(int index) { return _nshl_vec[index]; }

  /*! \brief Get the radius of the smallest sphere containing the particle.

   *

   * \param gc: `GeometryConfiguration *` Pointer to a `GeometryConfiguration` instance.

   * \return radius: `double` The radius of the sphere containing the particle.

   */

  double get_particle_radius(GeometryConfiguration *gc);

  /*! \brief Get the radius of a sphere by its index.

   *

   * This is a specialized function to get the radius of a sphere through its

  int proc_device;

  //! \brief Refinement mode selction flag.

  int refinemode;

  //! \brief flag defining a first iteration

  bool is_first_scale;

  //! \brief Maximum number of refinement iterations.

  int maxrefiters;

  //! \brief Required accuracy level.

  /*! \brief `ClusterIterationData` instance destroyer.

   */

  ~ClusterIterationData();

  /*! \brief Update field expansion orders.

   *

   * \param rcf: `double **` Matrix of sphere fractional radii.

   * \param inner_order: `int` The new inner expansion order to be set.

   * \param outer_order: `int` The new outer expansion order to be set.

   * \return result: `int` An exit code (0 if successful).

   */

  int update_orders(double** rcf, int inner_order, int outer_order);

};

// >>> END OF ClusterIterationData CLASS DEFINITION <<<

  int proc_device;

  //! \brief Refinement mode selction flag.

  int refinemode;

  //! \brief flag defining a first iteration

  bool is_first_scale;

  //! \brief Maximum number of refinement iterations.

  int maxrefiters;

  //! \brief Required accuracy level.

  /*! \brief `InclusionIterationData` instance destroyer.

   */

  ~InclusionIterationData();

  /*! \brief Update field expansion orders.

   *

   * \param rcf: `double **` Matrix of sphere fractional radii.

   * \param inner_order: `int` The new inner expansion order to be set.

   * \param outer_order: `int` The new outer expansion order to be set.

   * \return result: `int` An exit code (0 if successful).

   */

  int update_orders(double** rcf, int inner_order, int outer_order);

};

// >>> END OF InclusionIterationData CLASS DEFINITION <<< //

  double **tqss;

  //! \brief Scattering coefficients tensor.

  double ****zpv;

  //! \brief flag for first time initialisation

  bool is_first_scale;

  /*! \brief `SphereIterationData` default instance constructor.

   *

  /*! \brief `SphereIterationData` instance destroyer.

   */

  ~SphereIterationData();

  /*! \brief Update field expansion order.

   *

   * \param order: `int` The new expansion order to be set.

   * \return result: `int` An exit code (0 if successful).

   */

  int update_order(int order);

};

// >>> END OF SphereIterationData CLASS DEFINITION <<<

 * to radial coordinates, then it calls `sphar()` to calculate the vector of spherical

 * harmonics of the incident field.

 *

 * \param sconf: `ScattererConfiguration *` Pointer to scatterer configuration object.

 * \param rcf: `double **` Matrix of sphere configuration fractional radii.

 * \param c1: `ParticleDescriptor *` Pointer to a ParticleDescriptor instance.

 */

void str(ScattererConfiguration *sconf, ParticleDescriptor *c1);

void str(double **rcf, ParticleDescriptor *c1);

/*! \brief Compute radiation torques on particles on a k-vector oriented system.

 *