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

Implement dynamic order management in INCLUSION and fix wrong... 

Implement dynamic order management in INCLUSION and fix wrong ParticleDescriptorInclusion OpenMP copy constructor
parent a32b41ac

src/cluster/cluster.cpp

+1 −2
Original line number Diff line number Diff line
@@ -707,10 +707,9 @@ int cluster_jxi488_cycle( 
  const int npntts = gconf->npntts;

  const int isam = gconf->isam;

  const int jwtm = gconf->jwtm;

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

  int isq, ibf;

  int last_configuration;

  int num_configs = sconf->configurations;

  const int num_configs = sconf->configurations;

  const int ndirs = sa->nkks;

  int oindex = -1;

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

src/include/Commons.h

+2 −2
Original line number Diff line number Diff line
@@ -502,13 +502,13 @@ public: 
   */

  std::string get_descriptor_type() override { return "inclusion 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.

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

   */

  int update_orders(int inner_order, int outer_order) { return 0; }

  int update_orders(int inner_order, int outer_order);

};



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

src/include/IterationData.h

+9 −0
Original line number Diff line number Diff line
@@ -354,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 <<< //

src/include/inclu_subs.h

+2 −2
Original line number Diff line number Diff line
@@ -75,10 +75,10 @@ void indme( 


/*! \brief C++ porting of INSTR.

 *

 * \param sconf: `ScattererConfiguration *` Pointer to a ScattererConfiguration instance.

 * \param rcf: `double **` Pointer to matrix of fractional radii.

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

 */

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

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



/*! \brief C++ porting of OSPV.

 *

src/inclusion/inclusion.cpp

+100 −29
Original line number Diff line number Diff line
@@ -19,6 +19,7 @@ 
 * \brief Implementation of the calculation for a sphere with inclusions.

 */

#include <chrono>

#include <cmath>

#include <cstdio>

#include <exception>

#include <fstream>
@@ -233,12 +234,10 @@ void inclusion(const string& config_file, const string& data_file, const string& 
      // Open empty virtual ascii file for output

      InclusionOutputInfo *p_output = new InclusionOutputInfo(sconf, gconf, mpidata);

      InclusionIterationData *cid = new InclusionIterationData(gconf, sconf, mpidata, device_count);

      const np_int ndi = cid->c1->nsph * cid->c1->nlim;

      const np_int ndit = 2 * ndi;

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

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

      

      instr(sconf, cid->c1);

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

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

      double exdc = sconf->exdc;

      double exri = sqrt(exdc);
@@ -537,6 +536,7 @@ void inclusion(const string& config_file, const string& data_file, const string& 
    time_logger->log(message);

    fclose(timing_file);

    delete time_logger;

    delete logger;

  } // end instructions block of MPI process 0

  

    //===============================
@@ -673,8 +673,9 @@ void inclusion(const string& config_file, const string& data_file, const string& 
}



int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, InclusionIterationData *cid, InclusionOutputInfo *output, const string& output_path, VirtualBinaryFile *vtppoanp) {

  const dcomplex cc0 = 0.0 + I * 0.0;

  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);
@@ -684,24 +685,24 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
  int jer = 0;

  int lcalc = 0;

  int jaw = 1;

  int li = cid->c1->li;

  int le = cid->c1->le;

  int lm = cid->c1->lm;

  int nsph = cid->c1->nsph;

  bool is_first_scale = (jxi488 == 1);

  // int li = cid->c1->li;

  // int le = cid->c1->le;

  // int lm = cid->c1->lm;

  const int nsph = cid->c1->nsph;

  np_int mxndm = gconf->mxndm;

  int iavm = gconf->iavm;

  int inpol = gconf->in_pol;

  int npnt = cid->c1->npnt;

  int npntts = cid->c1->npntts;

  int isam = gconf->iavm;

  int jwtm = gconf->jwtm;

  np_int ndit = cid->c1->ndit;

  const int iavm = gconf->iavm;

  const int inpol = gconf->in_pol;

  const int npnt = cid->c1->npnt;

  const int npntts = cid->c1->npntts;

  const int isam = gconf->iavm;

  const int jwtm = gconf->jwtm;

  int isq, ibf;

  int last_configuration;

  dcomplex ent, entn;

  double enti;

  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;
@@ -725,6 +726,48 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
    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;

  // End of dynamic order check

  // label 120

  double sze = vkarg * cid->extr;

  last_configuration = 0;
@@ -745,7 +788,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
	  cid->c1->dc0[ic] = sconf->get_dielectric_constant(ic, i133 - 1, jxi488 - 1);

	// goes to 129

      } else { // label 127

	if (jxi488 == 1) {

	if (is_first_scale) {

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

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

	  }
@@ -921,7 +964,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
    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(
@@ -964,7 +1007,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
	    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)
@@ -974,7 +1017,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
		    );

	  }

	  // 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
@@ -1390,8 +1433,6 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
// >>> IMPLEMENTATION OF InclusionIterationData CLASS <<< //

InclusionIterationData::InclusionIterationData(GeometryConfiguration *gconf, ScattererConfiguration *sconf, const mixMPI *mpidata, const int device_count) {

  c1 = new ParticleDescriptorInclusion(gconf, sconf);

  const int ndi = c1->nsph * c1->nlim;

  const np_int ndit = 2 * ndi;

  gaps = new double[c1->nsph]();

  tqev = new double[3]();

  tqsv = new double[3]();
@@ -1502,8 +1543,6 @@ InclusionIterationData::InclusionIterationData(GeometryConfiguration *gconf, Sca 


InclusionIterationData::InclusionIterationData(const InclusionIterationData& rhs) {

  c1 = new ParticleDescriptorInclusion(reinterpret_cast<ParticleDescriptorInclusion &>(*(rhs.c1)));

  const int ndi = c1->nsph * c1->nlim;

  const np_int ndit = 2 * ndi;

  gaps = new double[c1->nsph]();

  for (int gi = 0; gi < c1->nsph; gi++) gaps[gi] = rhs.gaps[gi];

  tqev = new double[3]();
@@ -1654,8 +1693,6 @@ InclusionIterationData::InclusionIterationData(const InclusionIterationData& rhs 
#ifdef MPI_VERSION

InclusionIterationData::InclusionIterationData(const mixMPI *mpidata, const int device_count) {

  c1 = new ParticleDescriptorInclusion(mpidata);

  const int ndi = c1->nsph * c1->nlim;

  const np_int ndit = 2 * ndi;

  gaps = new double[c1->nsph]();

  MPI_Bcast(gaps, c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  tqev = new double[3]();
@@ -1792,8 +1829,6 @@ InclusionIterationData::InclusionIterationData(const mixMPI *mpidata, const int 


void InclusionIterationData::mpibcast(const mixMPI *mpidata) {

  c1->mpibcast(mpidata);

  const int ndi = c1->nsph * c1->nlim;

  const np_int ndit = 2 * ndi;

  MPI_Bcast(gaps, c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(tqev, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(tqsv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
@@ -1927,4 +1962,40 @@ InclusionIterationData::~InclusionIterationData() { 
  delete[] cmullr;

  delete[] cmul;

}



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

  int result = 0;

  int old_lm = c1->lm;

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

  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]();

      }

    }

  }

  instr(rcf, c1);

  thdps(c1->lm, zpv);

  delete[] am;

  delete[] am_vector;

  am_vector = new dcomplex[c1->ndm * c1->ndm]();

  am = new dcomplex*[c1->ndm];

  for (int ai = 0; ai < c1->ndm; ai++) {

    am[ai] = am_vector + ai * c1->ndm;

  }

  return result;

}

// >>> END OF InclusionIterationData CLASS IMPLEMENTATION <<<