Commit f2ec3daa authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Merge branch 'dynamic_orders' into 'master' 

Dynamic orders

See merge request giacomo.mulas/np_tmcode!94
parents a2047ebc f555f74f

src/cluster/cluster.cpp

+112 −25
Original line number Diff line number Diff line
@@ -19,6 +19,7 @@ 
 * \brief Implementation of the calculation for a cluster of spheres.

 */

#include <chrono>

#include <cmath>

#include <cstdio>

#include <exception>

#include <fstream>
@@ -228,6 +229,10 @@ void cluster(const string& config_file, const string& data_file, const string& o 
    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
@@ -240,7 +245,7 @@ void cluster(const string& config_file, const string& data_file, const string& o 
      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);
@@ -676,10 +681,14 @@ void cluster(const string& config_file, const string& data_file, const string& o 
#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);
@@ -689,24 +698,19 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
  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;

  bool is_first_scale = (jxi488 == 1);

  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;
@@ -716,7 +720,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
  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;
@@ -729,6 +733,49 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
    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;
@@ -754,7 +801,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	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);

	}
@@ -837,6 +884,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
    }

  }

#endif

  cid->refinemode = 0;

#ifdef DEBUG_AM

  VirtualAsciiFile *outam2 = new VirtualAsciiFile();

  string outam2_name = output_path + "/c_AM2_JXI" + to_string(jxi488) + ".txt";
@@ -972,7 +1020,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
    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(
@@ -984,7 +1032,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	  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);
@@ -1015,7 +1063,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	    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)
@@ -1025,7 +1073,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
		    );

	  }

	  // 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
@@ -2062,4 +2110,43 @@ ClusterIterationData::~ClusterIterationData() { 
  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 <<<

src/include/Commons.h

+23 −0
Original line number Diff line number Diff line
@@ -451,6 +451,14 @@ public: 
   * \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.
@@ -493,6 +501,14 @@ public: 
   * \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.
@@ -535,6 +551,13 @@ public: 
   * \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.

src/include/Configuration.h

+16 −3
Original line number Diff line number Diff line
@@ -288,11 +288,9 @@ protected: 
  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;
@@ -383,6 +381,8 @@ public: 
  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.

   *
@@ -511,6 +511,12 @@ public: 
   */

  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
@@ -521,6 +527,13 @@ public: 
   */

  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

src/include/IterationData.h

+25 −0
Original line number Diff line number Diff line
@@ -178,6 +178,15 @@ public: 
  /*! \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 <<<
@@ -345,6 +354,15 @@ public: 
  /*! \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 <<< //
@@ -503,6 +521,13 @@ public: 
  /*! \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 <<<

src/include/clu_subs.h

+2 −2
Original line number Diff line number Diff line
@@ -369,10 +369,10 @@ void scr2( 
 * 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.

 *