Commit 821dd790 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Move C2 data structures under ParticleDescriptor

parent 8ccbdb28

src/cluster/cluster.cpp

+6 −6
Original line number Diff line number Diff line
@@ -754,17 +754,17 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
      int ici = (nsh + 1) / 2;

      if (idfc == 0) {

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

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

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

      } else {

	if (jxi488 == 1) {

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

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

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

	}

      }

      if (nsh % 2 == 0) cid->c2->dc0[ici] = exdc;

      if (nsh % 2 == 0) cid->c1->dc0[ici] = exdc;

      dme(

	  cid->c1->li, i132, npnt, npntts, vkarg, exdc, exri,

	  cid->c1, cid->c2, jer, lcalc, cid->arg, last_configuration

	  cid->c1, jer, lcalc, cid->arg, last_configuration

	  );

      if (jer != 0) {

	sprintf(virtual_line, "  STOP IN DME\n");
@@ -915,10 +915,10 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
      sprintf(virtual_line, "     SPHERE %2d\n", i170);

      output->append_line(virtual_line);

      if (cid->c1->nshl[last_configuration - 1] != 1) {

	sprintf(virtual_line, "  SIZE=%15.7lE\n", cid->c2->vsz[i]);

	sprintf(virtual_line, "  SIZE=%15.7lE\n", cid->c1->vsz[i]);

	output->append_line(virtual_line);

      } else { // label 162

	sprintf(virtual_line, "  SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", cid->c2->vsz[i], real(cid->c2->vkt[i]), imag(cid->c2->vkt[i]));

	sprintf(virtual_line, "  SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", cid->c1->vsz[i], real(cid->c1->vkt[i]), imag(cid->c1->vkt[i]));

	output->append_line(virtual_line);

      }

      // label 164

src/include/Commons.h

+26 −58
Original line number Diff line number Diff line
@@ -40,62 +40,6 @@ 
class ParticleDescriptor;

class mixMPI;



/*! \brief Representation of the FORTRAN C2 blocks.

 *

 */

class C2 {

protected:

  //! \brief Number of spheres.

  int nsph;

  //! \brief Number of required orders.

  int nhspo;

  //! \brief QUESTION: what is nl?

  int nl;



public:

  //! \brief QUESTION: definition?

  dcomplex *ris;

  //! \brief QUESTION: definition?

  dcomplex *dlri;

  //! \brief Vector of dielectric constants.

  dcomplex *dc0;

  //! \brief QUESTION: definition?

  dcomplex *vkt;

  //! Vector of sphere sizes in units of 2*PI/LAMBDA

  double *vsz;



  /*! \brief C2 instance constructor.

   *

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

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

   */

  C2(GeometryConfiguration *gconf, ScattererConfiguration *sconf);



  /*! \brief C2 instance constructor copying its contents from preexisting instance.

   *

   * \param rhs: `C2` object to copy contents from

   */

  C2(const C2& rhs);



  //! \brief C2 instance destroyer.

  ~C2();



#ifdef MPI_VERSION

  /*! \brief C2 instance constructor copying all contents off MPI broadcast from MPI process 0

   *

   * \param mpidata: `mixMPI *` pointer to MPI data structure.

   */

  C2(const mixMPI *mpidata);



  /*! \brief send C2 instance from MPI process 0 via MPI broadcasts to all other processes

   *

   * \param mpidata: `mixMPI *` pointer to MPI data structure.

   */

  void mpibcast(const mixMPI *mpidata);

#endif



};



/*! \brief Representation of the FORTRAN C3 blocks.

 */

class C3 {
@@ -276,8 +220,6 @@ class ClusterIterationData { 
public:

  //! \brief Pointer to a C1 structure.

  ParticleDescriptor *c1;

  //! \brief Pointer to a C2 structure.

  C2 *c2;

  //! \brief Pointer to a C3 structure.

  C3 *c3;

  //! \brief Pointer to a C6 structure.
@@ -387,6 +329,14 @@ protected: 
  int _num_configurations;

  //! \brief Total number of layers from all sphere types.

  int _num_layers;

  //! \brief Space for different sphere types.

  int _nl;

  //! \brief NHSPO = 2 * MAX(NPNT,NPNTTS) - 1

  int _nhspo;

  //! \brief Number of points for numerical integration in layered spheres.

  int _npnt;

  //! \brief Number of points for numerical integration in transition layer.

  int _npntts;

  //! \brief Contiguous space for RMI.

  dcomplex *vec_rmi;

  //! \brief Contiguous space for REI.
@@ -462,6 +412,14 @@ public: 
  const int &num_configurations = _num_configurations;

  //! \brief Read-only view of total number of layers from all sphere types.

  const int &num_layers = _num_layers;

  //! \brief Read-only view on the space for different sphere configurations.

  const int &nl = _nl;

  //! \brief Read-only view of NHSPO.

  const int &nhspo = _nhspo;

  //! \brief Read-only view on number of points for numerical integration in layered spheres.

  const int &npnt = _npnt;

  //! \brief Read-only view on number of points for numerical integration in transition layer.

  const int &npntts = _npntts;

  //! \brief Matrix of inverse radial M coefficients.

  dcomplex **rmi;

  //! \brief Matrix of inverse radial E coefficients.
@@ -484,6 +442,16 @@ public: 
  int *iog;

  //! \brief Vector of number of layers in sphere type.

  int *nshl;

  //! \brief TBD

  dcomplex *ris;

  //! \brief TBD

  dcomplex *dlri;

  //! \brief TBD

  dcomplex *dc0;

  //! \brief TBD

  dcomplex *vkt;

  //! \brief Vector of sizes in units of 2*PI/LAMBDA

  double *vsz;

  // >>> END OF SECTION COMMON TO ALL DESCRIPTOR TYPES <<< //

  

  // >>> NEEDED BY SPHERE AND CLUSTER <<< //

src/include/inclu_subs.h

+1 −4
Original line number Diff line number Diff line
@@ -71,13 +71,10 @@ void incms(dcomplex **am, dcomplex **at, double enti, ParticleDescriptor *c1, C6 
 * \param lcalc: `int &` Maximum order achieved in calculation.

 * \param arg: `dcomplex &` TBD.

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

 * \param c2: `C2 *` Pointer to a C2 instance.

 */

void indme(

	   int i, int npnt, int npntts, double vk, dcomplex ent, double enti,

	   dcomplex entn, int &jer, int &lcalc, dcomplex &arg, ParticleDescriptor *c1,

	   C2 *c2

	   );

	   dcomplex entn, int &jer, int &lcalc, dcomplex &arg, ParticleDescriptor *c1	   );



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

 *

src/include/sph_subs.h

+4 −6
Original line number Diff line number Diff line
@@ -89,9 +89,8 @@ double cg1(int lmpml, int mu, int l, int m); 
 * \param ic: `int`

 * \param vk: `double`

 * \param c1: `ParticleDescriptor *` Pointer to `ParticleDescriptor` data structure.

 * \param c2: `C2 *` Pointer to `C2` data structure.

 */

void diel(int npntmo, int ns, int i, int ic, double vk, ParticleDescriptor *c1, C2 *c2);

void diel(int npntmo, int ns, int i, int ic, double vk, ParticleDescriptor *c1);



/*! \brief Compute Mie scattering coefficients.

 *
@@ -108,7 +107,6 @@ void diel(int npntmo, int ns, int i, int ic, double vk, ParticleDescriptor *c1, 
 * \param exdc: `double` External medium dielectric constant.

 * \param exri: `double` External medium refractive index.

 * \param c1: `ParticleDescriptor *` Pointer to a `ParticleDescriptor` data structure.

 * \param c2: `C2 *` Pointer to a `C2` data structure.

 * \param jer: `int &` Reference to integer error code variable.

 * \param lcalc: `int &` Reference to integer variable recording the maximum expansion order accounted for.

 * \param arg: `complex double &`
@@ -116,7 +114,7 @@ void diel(int npntmo, int ns, int i, int ic, double vk, ParticleDescriptor *c1, 
 */

void dme(

	 int li, int i, int npnt, int npntts, double vk, double exdc, double exri,

	 ParticleDescriptor *c1, C2 *c2, int &jer, int &lcalc, dcomplex &arg, int last_conf=0

	 ParticleDescriptor *c1, int &jer, int &lcalc, dcomplex &arg, int last_conf=0

);



/*! \brief Bessel function calculation control parameters.
@@ -262,11 +260,11 @@ void rkc( 
 * \param y2: `complex double &`

 * \param dy1: `complex double &`

 * \param dy2: `complex double &`

 * \param c2: `C2 *` Pointer to a `C2` data structure.

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

 */

void rkt(

	 int npntmo, double step, double &x, int lpo, dcomplex &y1,

	 dcomplex &y2, dcomplex &dy1, dcomplex &dy2, C2 *c2

	 dcomplex &y2, dcomplex &dy1, dcomplex &dy2, ParticleDescriptor *c1

);



/*! \brief Spherical Bessel functions.

src/libnptm/Commons.cpp

+60 −77
Original line number Diff line number Diff line
@@ -36,78 +36,6 @@ 
#include <mpi.h>

#endif



C2::C2(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {

  nsph = gconf->number_of_spheres;

  int npnt = gconf->npnt;

  int npntts = gconf->npntts;

  int max_n = (npnt > npntts) ? npnt : npntts;

  nhspo = 2 * max_n - 1;

  nl = sconf->configurations;

  if (nsph == 1 && nl == 1) nl = 5;

  ris = new dcomplex[nhspo]();

  dlri = new dcomplex[nhspo]();

  vkt = new dcomplex[nsph]();

  dc0 = new dcomplex[nl]();

  vsz = new double[nsph]();

}



C2::C2(const C2& rhs) {

  nsph = rhs.nsph;

  nhspo = rhs.nhspo;

  nl = rhs.nl;

  ris = new dcomplex[nhspo]();

  dlri = new dcomplex[nhspo]();

  for (int ind=0; ind<nhspo; ind++) {

    ris[ind] = rhs.ris[ind];

    dlri[ind] = rhs.dlri[ind];

  }

  vkt = new dcomplex[nsph]();

  vsz = new double[nsph]();

  for (int ind=0; ind<nsph; ind++) {

    vkt[ind] = rhs.vkt[ind];

    vsz[ind] = rhs.vsz[ind];

  }

  dc0 = new dcomplex[nl]();

  for (int ind=0; ind<nl; ind++) dc0[ind] = rhs.dc0[ind];

}



C2::~C2() {

  delete[] ris;

  delete[] dlri;

  delete[] vkt;

  delete[] dc0;

  delete[] vsz;

}



#ifdef MPI_VERSION

C2::C2(const mixMPI *mpidata) {

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

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

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

  ris = new dcomplex[nhspo]();

  dlri = new dcomplex[nhspo]();

  MPI_Bcast(ris, nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(dlri, nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  vkt = new dcomplex[nsph]();

  vsz = new double[nsph]();

  MPI_Bcast(vkt, nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(vsz, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  dc0 = new dcomplex[nl]();

  MPI_Bcast(dc0, nl, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

}



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

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

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

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

  MPI_Bcast(ris, nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(dlri, nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(vkt, nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(vsz, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(dc0, nl, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

}

#endif



C3::C3() {

  tsas = new dcomplex*[2];

  tsas[0] = new dcomplex[2];
@@ -306,7 +234,6 @@ mixMPI::~mixMPI() { 


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

  c1 = new ParticleDescriptorCluster(gconf, sconf);

  c2 = new C2(gconf, sconf);

  c3 = new C3();

  c6 = new C6(c1->lmtpo);

  const int ndi = c1->nsph * c1->nlim;
@@ -411,7 +338,6 @@ ClusterIterationData::ClusterIterationData(GeometryConfiguration *gconf, Scatter 


ClusterIterationData::ClusterIterationData(const ClusterIterationData& rhs) {

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

  c2 = new C2(*(rhs.c2));

  c3 = new C3(*(rhs.c3));

  c6 = new C6(*(rhs.c6));

  const int ndi = c1->nsph * c1->nlim;
@@ -560,7 +486,6 @@ ClusterIterationData::ClusterIterationData(const ClusterIterationData& rhs) { 
#ifdef MPI_VERSION

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

  c1 = new ParticleDescriptorCluster(mpidata);

  c2 = new C2(mpidata);

  c3 = new C3(mpidata);

  c6 = new C6(mpidata);

  const int ndi = c1->nsph * c1->nlim;
@@ -697,7 +622,6 @@ ClusterIterationData::ClusterIterationData(const mixMPI *mpidata, const int devi 


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

  c1->mpibcast(mpidata);

  c2->mpibcast(mpidata);

  c3->mpibcast(mpidata);

  c6->mpibcast(mpidata);

  const int ndi = c1->nsph * c1->nlim;
@@ -783,7 +707,6 @@ ClusterIterationData::~ClusterIterationData() { 
  }

  delete[] zpv;

  delete c1;

  delete c2;

  delete c3;

  delete c6;

  delete c9;
@@ -906,6 +829,18 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo 
    ros[ci] = sconf->get_radius(ci);

    nshl[ci] = sconf->get_nshl(ci);

  }

  _npnt = gconf->npnt;

  _npntts = gconf->npntts;

  int max_n = (npnt > npntts) ? npnt : npntts;

  _nhspo = 2 * max_n - 1;

  _nl = sconf->configurations;

  if (_nsph == 1 && _nl == 1) _nl = 5;

  ris = new dcomplex[_nhspo]();

  dlri = new dcomplex[_nhspo]();

  vkt = new dcomplex[_nsph]();

  dc0 = new dcomplex[_nl]();

  vsz = new double[_nsph]();

  

  // >>> NEEDED BY SPHERE AND CLUSTER <<<

  sas = NULL;

  vints = NULL;
@@ -1021,6 +956,26 @@ ParticleDescriptor::ParticleDescriptor(const ParticleDescriptor &rhs) { 
    ros[ci] = rhs.ros[ci];

    nshl[ci] = rhs.nshl[ci];

  }

  _npnt = rhs._npnt;

  _npntts = rhs._npntts;

  _nhspo = rhs._nhspo;

  _nl = rhs._nl;

  ris = new dcomplex[_nhspo]();

  dlri = new dcomplex[_nhspo]();

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

    ris[ri] = rhs.ris[ri];

    dlri[ri] = rhs.dlri[ri];

  }

  vkt = new dcomplex[_nsph]();

  vsz = new double[_nsph]();

  for (int vi = 0; vi < _nsph; vi++) {

    vkt[vi] = rhs.vkt[vi];

    vsz[vi] = rhs.vsz[vi];

  }

  dc0 = new dcomplex[_nl]();

  for (int di = 0; di < _nl; di++) {

    dc0[di] = rhs.dc0[di];

  }

  // >>> NEEDED BY SPHERE AND CLUSTER <<<

  sas = NULL;

  vints = NULL;
@@ -1131,6 +1086,20 @@ ParticleDescriptor::ParticleDescriptor(const mixMPI *mpidata) { 
  MPI_Bcast(iog, _nsph, MPI_INT, 0, MPI_COMM_WORLD);

  ros = new double[_num_configurations];

  MPI_Bcast(ros, _num_configurations, MPI_DOUBLE, 0, MPI_COMM_WORLD);

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

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

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

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

  ris = new dcomplex[_nhspo];

  MPI_Bcast(ris, _nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  dlri = new dcomplex[_nhspo];

  MPI_Bcast(dlri, _nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  vkt = new dcomplex[_nsph];

  MPI_Bcast(vkt, _nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  vsz = new double[_nsph];

  MPI_Bcast(vsz, _nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  dc0 = new dcomplex[_nl];

  MPI_Bcast(dc0, _nl, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

}



void ParticleDescriptor::mpibcast(const mixMPI *mpidata) {
@@ -1163,6 +1132,15 @@ void ParticleDescriptor::mpibcast(const mixMPI *mpidata) { 
  MPI_Bcast(rzz, _nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(iog, _nsph, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(ros, _num_configurations, MPI_DOUBLE, 0, MPI_COMM_WORLD);

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

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

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

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

  MPI_Bcast(ris, _nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(dlri, _nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(vkt, _nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(vsz, _nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(dc0, _nl, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  // >>> NEEDED BY SPHERE AND CLUSTER <<< //

  if (_class_type == SPHERE_TYPE || _class_type == CLUSTER_TYPE) {

    MPI_Bcast(vec_sas, 4 * _nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
@@ -1237,6 +1215,11 @@ ParticleDescriptor::~ParticleDescriptor() { 
  delete[] vec_rei;

  delete[] rmi;

  delete[] vec_rmi;

  delete[] ris;

  delete[] dlri;

  delete[] vkt;

  delete[] dc0;

  delete[] vsz;

  // Inclusion class members, destroyed only if sub-class is INCLUSION

  if (_class_type == INCLUSION_TYPE) {

    delete[] rm0;