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

Implement dynamic order management in SPHERE

parent 2f9ec364

src/include/Commons.h

+4 −4
Original line number Diff line number Diff line
@@ -452,7 +452,7 @@ public: 
   */

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

  

  /*! \brief Interface function to update field expansion orders.

  /*! \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.
@@ -552,12 +552,12 @@ public: 
   */

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

  

  /*! \brief Interface function to update field expansion orders.

  /*! \brief Update the field expansion order.

   *

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

   * \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) { return 0; }

  int update_order(int order);

};



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

src/include/IterationData.h

+7 −0
Original line number Diff line number Diff line
@@ -512,6 +512,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/libnptm/Commons.cpp

+27 −0
Original line number Diff line number Diff line
@@ -1397,6 +1397,33 @@ ParticleDescriptorSphere::ParticleDescriptorSphere(const mixMPI *mpidata) : Part 
  }

}

#endif // MPI_VERSION



int ParticleDescriptorSphere::update_order(int order) {

  int result = 0;

  if (order != _lm) {

    _lm = order;

    _li = order;

    delete[] rmi;

    delete[] vec_rmi;

    delete[] rei;

    delete[] vec_rei;

    vec_rmi = new dcomplex[_li * _nsph]();

    rmi = new dcomplex*[_li];

    vec_rei = new dcomplex[_li * _nsph]();

    rei = new dcomplex*[_li];

    for (int ri = 0; ri < _li; ri++) {

      rmi[ri] = vec_rmi + (_nsph * ri);

      rei[ri] = vec_rei + (_nsph * ri);

    }

    int nlmmt = 2 * _nsph * _li * (_li + 2);

    delete[] w;

    delete[] vec_w;

    vec_w = new dcomplex[nlmmt * 4]();

    w = new dcomplex*[nlmmt];

    for (int wi = 0; wi < nlmmt; wi++) w[wi] = vec_w + (4 * wi);

  }

  return result;

}

// >>> End of ParticleDescriptorSphere class implementation. <<< //



// >>> ScatteringAngles class implementation. <<< //

src/sphere/sphere.cpp

+59 −9
Original line number Diff line number Diff line
@@ -542,11 +542,13 @@ int sphere_jxi488_cycle( 
  const double pi = 2.0 * half_pi;

  int jer = 0;

  Logger *logger = new Logger(LOG_INFO);

  string message = "INIT";

  int oindex = 0;

  int jxi = jxi488 + 1;

  int jxindex = jxi - oi->first_xi;

  bool is_first_scale = (jxi == 1);

  int nsph = gconf->number_of_spheres;

  int l_max = gconf->l_max;

  //int l_max = gconf->l_max;

  int in_pol = gconf->in_pol;

  int npnt = gconf->npnt;

  int npntts = gconf->npntts;
@@ -576,12 +578,30 @@ int sphere_jxi488_cycle( 
    oi->vec_vk[jxindex] = vk;

    oi->vec_xi[jxindex] = xi;

  }

  // Adaptive definition of L_MAX

  double wavelength = 2.0 * pi / vk;

  double size_param = 2.0 * pi * sconf->get_radius(0) / wavelength;

  int N = int(size_param + 4.05 * pow(size_param, 1.0 / 3.0)) + 2;

  if (N < l_max) l_max = N;

  // End of adaptive definition of L_MAX

  // Dynamic order check

  const int max_lm = gconf->l_max;

  int l_max = gconf->l_max;

  const double alamb = 2.0 * pi / vk;

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

  int recommended_lm = 4 + (int)ceil(size_par_lm + 4.05 * pow(size_par_lm, 1.0 / 3.0));

  if (recommended_lm != l_max) {

    if (recommended_lm > max_lm) {

      message = "WARNING: internal order " + to_string(max_lm) + " for scale iteration "

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

	+ ").\n";

      logger->log(message, LOG_WARN);

    } else if (recommended_lm < max_lm) {

      int new_lm = recommended_lm;

      message = "INFO: lowering internal order from " + to_string(max_lm) + " to "

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

      logger->log(message, LOG_INFO);

      sid->update_order(new_lm);

      is_first_scale = true;

      jw = 1;

      l_max = new_lm;

    }

  }

  // End of dynamic order check

  vtppoanp->append_line(VirtualBinaryLine(vk));

  double thsca = (gconf->isam > 1) ? sa->ths - sa->th : 0.0;

  for (int i132 = 0; i132 < nsph; i132++) {
@@ -601,7 +621,7 @@ int sphere_jxi488_cycle( 
      for (int ic = 0; ic < ici; ic++)

	sid->c1->dc0[ic] = sconf->get_dielectric_constant(ic, i132, jxi488); // WARNING: IDFC=0 is not tested!

    } else { // IDFC != 0

      if (jxi == 1) {

      if (is_first_scale) {

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

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

	}
@@ -766,7 +786,7 @@ int sphere_jxi488_cycle( 
	      wmamp(2, sid->costs, sid->sints, sid->cosps, sid->sinps, in_pol, l_max, 0, nsph, sid->args, sid->us, sid->upsmp, sid->unsmp, sid->c1);

	    }

	  }

	  if (nkks != 0 || jxi == 1) {

	  if (nkks != 0 || is_first_scale) {

	    upvsp(

	      sid->u, sid->upmp, sid->unmp, sid->us, sid->upsmp, sid->unsmp,

	      sid->up, sid->un, sid->ups, sid->uns, sid->duk, isq, ibf,
@@ -1253,4 +1273,34 @@ int SphereIterationData::mpibcast(const mixMPI *mpidata) { 
  return 0;

}

#endif // MPI_VERSION



int SphereIterationData::update_order(int order) {

  int result = 0;

  int old_lm = _lm;

  if (order != _lm) {

    _lm = order;

    ((ParticleDescriptorSphere *)c1)->update_order(order);

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

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

	delete[] zpv[zi][zj];

      }

      delete[] zpv[zi];

    }

    delete[] zpv;

    delete[] vec_zpv;

    vec_zpv = new double[_lm * 12]();

    zpv = new double***[_lm];

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

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

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

	int vec_index = 12 * zi + 4 * zj;

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

	zpv[zi][zj][0] = (vec_zpv + vec_index);

	zpv[zi][zj][1] = (vec_zpv + vec_index + 2);

      }

    }

    thdps(_lm, zpv);

  }

  return result;

}

// >>> END OF SphereIterationData CLASS IMPLEMENTATION <<<