diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index f3e9141dd5567e96a0155c03099870a1f6d7fcd8..f9e43625f7c346bf32dedc9661d734002346f65b 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -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);
@@ -283,53 +288,55 @@ void cluster(const string& config_file, const string& data_file, const string& o
}
// 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;
@@ -420,7 +427,8 @@ void cluster(const string& config_file, const string& data_file, const string& o
}
#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;
@@ -468,6 +476,10 @@ void cluster(const string& config_file, const string& data_file, const string& o
// 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];
@@ -599,7 +611,8 @@ void cluster(const string& config_file, const string& data_file, const string& o
// 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;
@@ -676,10 +689,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 +706,20 @@ 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;
+ // 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;
@@ -716,7 +729,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 +742,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 +810,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 +893,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 +1029,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 +1041,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 +1072,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 +1082,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
@@ -1515,6 +1572,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
} // 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
@@ -1631,6 +1690,8 @@ ClusterIterationData::ClusterIterationData(GeometryConfiguration *gconf, Scatter
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
@@ -1779,6 +1840,7 @@ ClusterIterationData::ClusterIterationData(const ClusterIterationData& rhs) {
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;
@@ -1916,6 +1978,7 @@ ClusterIterationData::ClusterIterationData(const mixMPI *mpidata, const int devi
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);
}
@@ -1987,6 +2050,7 @@ void ClusterIterationData::mpibcast(const mixMPI *mpidata) {
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);
}
@@ -2062,4 +2126,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 <<<
diff --git a/src/include/Commons.h b/src/include/Commons.h
index 31efcc39feb264b0b62b275076477b9402a5e475..54f062e5d35adc5a808193ed6b1441d23ecac927 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -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.
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
index d23aa6ff178ee60c8f251ab9b88a86349b6e1762..4235f56d9d42cb4849380282f66b90871b5bf1cb 100644
--- a/src/include/Configuration.h
+++ b/src/include/Configuration.h
@@ -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
diff --git a/src/include/IterationData.h b/src/include/IterationData.h
index bf38b033096cfd0329dd53ddf075a2f70f36ec26..a8d82dbf9e2d5a02d583fc94cf3f5f80442437dd 100644
--- a/src/include/IterationData.h
+++ b/src/include/IterationData.h
@@ -135,6 +135,8 @@ public:
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.
@@ -178,6 +180,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 <<<
@@ -302,6 +313,8 @@ public:
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.
@@ -345,6 +358,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 <<< //
@@ -464,6 +486,8 @@ public:
double **tqss;
//! \brief Scattering coefficients tensor.
double ****zpv;
+ //! \brief flag for first time initialisation
+ bool is_first_scale;
/*! \brief `SphereIterationData` default instance constructor.
*
@@ -503,6 +527,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 <<<
diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index 9c5feffe260086c5b8a42fa7608a5bb82585f35a..40c5f0e322702b1a3e3c7801cdaf400fc792e852 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -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.
*
diff --git a/src/include/inclu_subs.h b/src/include/inclu_subs.h
index 64e5a1f3da5d11fc060f709b60ac50d4d2055554..c7af5e26c70594e4a0d9184acf33fffc1de9113e 100644
--- a/src/include/inclu_subs.h
+++ b/src/include/inclu_subs.h
@@ -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.
*
diff --git a/src/inclusion/inclusion.cpp b/src/inclusion/inclusion.cpp
index 0b0033c4a20f7609c5fd5e9dc940efdf8f7509df..b70306c044b424f3d719d8df501b98f605964d44 100644
--- a/src/inclusion/inclusion.cpp
+++ b/src/inclusion/inclusion.cpp
@@ -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);
@@ -279,50 +278,50 @@ void inclusion(const string& config_file, const string& data_file, const string&
}
// do the first iteration on jxi488 separately, since it seems to be different from the others
- int jxi488 = 1;
+ 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 = inclusion_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, cid, p_output, output_path, vtppoanp);
- }
- }
-#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 = inclusion_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, cid, p_output, output_path, vtppoanp);
+// }
+// }
+// #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;
- /* 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;
- }
+ // 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;
+ // }
//==================================================
// do the first outputs here, so that I open here the new files, afterwards I only append
@@ -416,7 +415,7 @@ void inclusion(const string& config_file, const string& data_file, const string&
}
#pragma omp barrier
// ok, now I can actually start the parallel calculations
- for (int ixi488=2; ixi488<=cid_2->number_of_scales; ixi488 +=myMPIstride) {
+ 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;
@@ -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
//===============================
@@ -595,7 +595,7 @@ void inclusion(const string& config_file, const string& data_file, const string&
// 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) {
+ 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;
@@ -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,26 @@ 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;
+ // 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;
+ //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 +728,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 +790,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 +966,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 +1009,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 +1019,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
@@ -1371,6 +1416,8 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
} // 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
@@ -1390,8 +1437,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]();
@@ -1493,6 +1538,8 @@ InclusionIterationData::InclusionIterationData(GeometryConfiguration *gconf, Sca
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
@@ -1502,8 +1549,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]();
@@ -1646,6 +1691,7 @@ InclusionIterationData::InclusionIterationData(const InclusionIterationData& rhs
extr = rhs.extr;
proc_device = rhs.proc_device;
+ is_first_scale = rhs.is_first_scale;
refinemode = rhs.refinemode;
maxrefiters = rhs.maxrefiters;
accuracygoal = rhs.accuracygoal;
@@ -1654,8 +1700,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]();
@@ -1786,14 +1830,13 @@ InclusionIterationData::InclusionIterationData(const mixMPI *mpidata, const int
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);
}
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);
@@ -1853,6 +1896,7 @@ void InclusionIterationData::mpibcast(const mixMPI *mpidata) {
MPI_Bcast(&nimd, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(&extr, 1, MPI_DOUBLE, 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);
}
@@ -1927,4 +1971,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 <<<
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index 29bf18b637b149e7d3d6b9fa71ae8433c7e135cd..b93abab2ce2fa19c341256ed4a9847901eb99f4c 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -83,10 +83,12 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo
_ndm = 0;
gcs = 0.0;
_num_configurations = sconf->configurations;
- _num_layers = (sconf->use_external_sphere) ? 1 : 0;
+ // _num_layers = (sconf->use_external_sphere) ? 1 : 0;
+ _num_layers = 0;
_max_layers = 1;
for (int nli = 0; nli < num_configurations; nli++) {
int nl = sconf->get_nshl(nli);
+ if (nli == 0 && sconf->use_external_sphere) nl++;
_num_layers += nl;
if (nl > _max_layers) _max_layers = nl;
}
@@ -136,8 +138,6 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo
_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]();
@@ -310,7 +310,6 @@ ParticleDescriptor::ParticleDescriptor(const ParticleDescriptor &rhs) {
tfsas = 0.0 + I * 0.0;
vec_tsas = NULL;
tsas = NULL;
- gcs = 0.0;
scs = 0.0;
ecs = 0.0;
acs = 0.0;
@@ -982,6 +981,91 @@ ParticleDescriptorCluster::ParticleDescriptorCluster(const mixMPI *mpidata) : Pa
MPI_Bcast(rac3j, _lmtpo, MPI_DOUBLE, 0, MPI_COMM_WORLD);
}
#endif // MPI_VERSION
+
+int ParticleDescriptorCluster::update_orders(int inner_order, int outer_order) {
+ int result = 0;
+ bool changed_li = false;
+ bool changed_le = false;
+ if (inner_order != _li) {
+ _li = inner_order;
+ changed_li = true;
+ delete[] vec_rmi;
+ delete[] rmi;
+ vec_rmi = new dcomplex[_li * _nsph]();
+ rmi = new dcomplex*[_li];
+ delete[] vec_rei;
+ delete[] rei;
+ 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);
+ }
+ _litpo = _li + _li + 1;
+ _litpos = _litpo * _litpo;
+ delete[] vh;
+ vh = new dcomplex[_ncou * _litpo]();
+ delete[] vyhj;
+ vyhj = new dcomplex[_ncou * _litpos]();
+ _nlim = _li * (_li + 2);
+ _ndi = _nsph * _nlim;
+ _ndit = 2 * _nsph * _nlim;
+ }
+ if (outer_order != _le) {
+ _le = outer_order;
+ changed_le = true;
+ _nlem = _le * (_le + 2);
+ _nlemt = 2 * _nlem;
+ delete[] vec_am0m;
+ vec_am0m = new dcomplex[_nlemt * _nlemt]();
+ delete[] am0m;
+ am0m = new dcomplex*[_nlemt];
+ for (int ai = 0; ai < _nlemt; ai++) am0m[ai] = vec_am0m + (ai * _nlemt);
+ delete[] vec_w;
+ delete[] w;
+ vec_w = new dcomplex[_nlemt * 4]();
+ w = new dcomplex*[nlemt];
+ for (int wi = 0; wi < nlemt; wi++) w[wi] = vec_w + (4 * wi);
+ }
+ if (changed_li || changed_le) {
+ _lm = (_li > _le) ? _li : _le;
+ _lmpo = _lm + 1;
+ _lmtpo = _li + _le + 1;
+ _lmtpos = _lmtpo * _lmtpo;
+ _nv3j = (_lm * (_lm + 1) * (2 * _lm + 7)) / 6;
+ delete[] vec_ind3j;
+ vec_ind3j = new int[(_lm + 1) * _lm]();
+ delete[] ind3j;
+ ind3j = new int*[_lm + 1];
+ for (int ii = 0; ii <= _lm; ii++) ind3j[ii] = vec_ind3j + (_lm * ii);
+ delete[] vj0;
+ vj0 = new dcomplex[_nsph * _lmtpo]();
+ delete[] vyj0;
+ vyj0 = new dcomplex[_nsph * _lmtpos]();
+ delete[] v3j0;
+ v3j0 = new double[_nv3j]();
+ delete[] rac3j;
+ rac3j = new double[_lmtpo]();
+ delete[] vec_gis;
+ vec_gis = new dcomplex[_ndi * _nlem]();
+ delete[] gis;
+ gis = new dcomplex*[_ndi];
+ delete[] vec_gls;
+ vec_gls = new dcomplex[_ndi * _nlem]();
+ delete[] gls;
+ gls = new dcomplex*[_ndi];
+ for (int gi = 0; gi < _ndi; gi++) {
+ gis[gi] = vec_gis + (gi * _nlem);
+ gls[gi] = vec_gls + (gi * _nlem);
+ }
+ delete[] vec_sam;
+ vec_sam = new dcomplex[_ndit * _nlemt]();
+ delete[] sam;
+ sam = new dcomplex*[_ndit];
+ for (int si = 0; si < _ndit; si++) sam[si] = vec_sam + (si * _nlemt);
+ }
+ return result;
+}
// >>> End of ParticleDescriptorCluster class implementation. <<< //
// >>> ParticleDescriptorInclusion class implementation. <<< //
@@ -1066,6 +1150,8 @@ ParticleDescriptorInclusion::ParticleDescriptorInclusion(const ParticleDescripto
vec_ind3j = new int[(_lm + 1) * _lm];
np_int vec_ind3j_size = (_lm + 1) * _lm;
for (np_int ii = 0; ii < vec_ind3j_size; ii++) vec_ind3j[ii] = rhs.vec_ind3j[ii];
+ ind3j = new int*[_lm + 1];
+ for (int ii = 0; ii <= _lm; ii++) ind3j[ii] = vec_ind3j + (_lm * ii);
vh = new dcomplex[_ncou * _litpo];
for (int hi = 0; hi < _ncou * _litpo; hi++) vh[hi] = rhs.vh[hi];
@@ -1106,9 +1192,7 @@ ParticleDescriptorInclusion::ParticleDescriptorInclusion(const ParticleDescripto
ecscm[ci] = rhs.ecscm[ci];
}
v3j0 = new double[_nv3j];
- for (int vj = 0; vj < _nv3j; vj++) v3j0[vj] = rhs.v3j0[_nv3j];
- ind3j = new int*[_lm + 1];
- for (int ii = 0; ii <= _lm; ii++) ind3j[ii] = vec_ind3j + (_lm * ii);
+ for (int vj = 0; vj < _nv3j; vj++) v3j0[vj] = rhs.v3j0[vj];
rac3j = new double[_lmtpo];
for (int ri = 0; ri < _lmtpo; ri++) rac3j[ri] = rhs.rac3j[ri];
// >>> NEEDED BY INCLU <<< //
@@ -1217,6 +1301,96 @@ ParticleDescriptorInclusion::ParticleDescriptorInclusion(const mixMPI *mpidata)
for (int ai = 0; ai < _nlemt; ai++) at[ai] = vec_at + (ai * _ndm);
}
#endif // MPI_VERSION
+
+int ParticleDescriptorInclusion::update_orders(int inner_order, int outer_order) {
+ int result = 0;
+ bool changed_li = false;
+ bool changed_le = false;
+ if (inner_order != _li) {
+ _li = inner_order;
+ changed_li = true;
+ _nlim = _li * (_li + 2);
+ _litpo = _li + _li + 1;
+ _litpos = _litpo * _litpo;
+ _ndi = _nsph * _nlim;
+ _ndit = 2 * _nsph * _nlim;
+ 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);
+ }
+ delete[] vh;
+ vh = new dcomplex[_ncou * _litpo]();
+ delete[] vyhj;
+ vyhj = new dcomplex[_ncou * _litpos]();
+ }
+ if (outer_order != le) {
+ _le = outer_order;
+ changed_le = true;
+ _nlem = _le * (_le + 2);
+ _nlemt = 2 * _nlem;
+ delete[] w;
+ delete[] vec_w;
+ vec_w = new dcomplex[_nlemt * 4]();
+ w = new dcomplex*[_nlemt];
+ for (int wi = 0; wi < _nlemt; wi++) w[wi] = vec_w + (4 * wi);
+ delete[] am0m;
+ delete[] vec_am0m;
+ vec_am0m = new dcomplex[_nlemt * _nlemt]();
+ am0m = new dcomplex*[_nlemt];
+ for (int ai = 0; ai < _nlemt; ai++) am0m[ai] = vec_am0m + (ai * _nlemt);
+ delete[] rm0;
+ rm0 = new dcomplex[_le]();
+ delete[] re0;
+ re0 = new dcomplex[_le]();
+ delete[] rmw;
+ rmw = new dcomplex[_le]();
+ delete[] rew;
+ rew = new dcomplex[_le]();
+ delete[] tm;
+ tm = new dcomplex[_le]();
+ delete[] te;
+ te = new dcomplex[_le]();
+ delete[] tm0;
+ tm0 = new dcomplex[_le]();
+ delete[] te0;
+ te0 = new dcomplex[_le]();
+ }
+ if (changed_li || changed_le) {
+ _lm = (_li > _le) ? _li : _le;
+ _lmpo = _lm + 1;
+ _lmtpo = _li + _le + 1;
+ _lmtpos = _lmtpo * _lmtpo;
+ _nv3j = (_lm * (_lm + 1) * (2 * _lm + 7)) / 6;
+ _ndm = 2 * (_nsph * _nlim + _nlem);
+ delete[] ind3j;
+ delete[] vec_ind3j;
+ vec_ind3j = new int[(_lm + 1) * _lm]();
+ ind3j = new int*[_lm + 1];
+ for (int ii = 0; ii <= _lm; ii++) ind3j[ii] = vec_ind3j + (_lm * ii);
+ delete[] vj0;
+ vj0 = new dcomplex[_nsph * _lmtpo]();
+ delete[] vyj0;
+ vyj0 = new dcomplex[_nsph * _lmtpos]();
+ delete[] v3j0;
+ v3j0 = new double[_nv3j]();
+ delete[] rac3j;
+ rac3j = new double[_lmtpo]();
+ delete[] at;
+ delete[] vec_at;
+ vec_at = new dcomplex[_nlemt * _ndm]();
+ at = new dcomplex*[_nlemt];
+ for (int ai = 0; ai < _nlemt; ai++) at[ai] = vec_at + (ai * _ndm);
+ }
+ return result;
+}
// >>> End of ParticleDescriptorInclusion class implementation. <<< //
// >>> ParticleDescriptorSphere class implementation. <<< //
@@ -1312,6 +1486,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. <<< //
diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp
index 68e076630b26c5eb8bb91b2a1dcd5d8b392392d5..2aa61343aaf91e99294d2d786c72ddae2f290c5b 100644
--- a/src/libnptm/Configuration.cpp
+++ b/src/libnptm/Configuration.cpp
@@ -971,6 +971,40 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(const std::string& f
return conf;
}
+double ScattererConfiguration::get_max_radius() {
+ double result = 0.0;
+ for (int ci = 0; ci < _configurations; ci++) {
+ if (_radii_of_spheres[ci] > result)
+ result = _radii_of_spheres[ci];
+ }
+ return result;
+}
+
+double ScattererConfiguration::get_particle_radius(GeometryConfiguration *gc) {
+ double result = 0.0;
+ if (_use_external_sphere) {
+ result = _radii_of_spheres[0] * _rcf[0][_nshl_vec[0] - 1];
+ } else {
+ double x, y, z;
+ int far_index = -1;
+ double dist2, max_dist;
+ double max_dist2 = 0.0;
+ for (int si = 0; si < _number_of_spheres; si++) {
+ x = gc->get_sph_x(si);
+ y = gc->get_sph_y(si);
+ z = gc->get_sph_z(si);
+ dist2 = x * x + y * y + z * z;
+ if (dist2 > max_dist2) {
+ max_dist2 = dist2;
+ far_index = si;
+ }
+ }
+ max_dist = sqrt(max_dist2);
+ result = max_dist + _radii_of_spheres[far_index];
+ }
+ return result;
+}
+
void ScattererConfiguration::print() {
int ies = (_use_external_sphere)? 1 : 0;
printf("### CONFIGURATION DATA ###\n");
diff --git a/src/libnptm/clu_subs.cpp b/src/libnptm/clu_subs.cpp
index 0199f31a96dd9d759be40f90c07d70cef76c7e0b..3548af3d656cff48a18fddafafb2e61a04542852 100644
--- a/src/libnptm/clu_subs.cpp
+++ b/src/libnptm/clu_subs.cpp
@@ -2442,7 +2442,7 @@ void scr2(
#endif
}
-void str(ScattererConfiguration *sconf, ParticleDescriptor *c1) {
+void str(double **rcf, ParticleDescriptor *c1) {
dcomplex *ylm;
const double pi = acos(-1.0);
int last_configuration;
@@ -2457,7 +2457,7 @@ void str(ScattererConfiguration *sconf, ParticleDescriptor *c1) {
c1->gcsv[i18 - 1] = gcss;
int nsh = c1->nshl[last_configuration - 1];
for (int j16 = 1; j16 <= nsh; j16++) {
- c1->rc[last_configuration - 1][j16 - 1] = sconf->get_rcf(last_configuration - 1, j16 - 1) * c1->ros[last_configuration - 1];
+ c1->rc[last_configuration - 1][j16 - 1] = rcf[last_configuration - 1][j16 - 1] * c1->ros[last_configuration - 1];
} // j16 loop
}
c1->gcs += gcss;
diff --git a/src/libnptm/inclu_subs.cpp b/src/libnptm/inclu_subs.cpp
index e2eede6927f1701ca34656ba8cd1a6051baeecdc..12859ab2a1953eb2387fa5d3c9eb804fbb8c4f0b 100644
--- a/src/libnptm/inclu_subs.cpp
+++ b/src/libnptm/inclu_subs.cpp
@@ -493,7 +493,7 @@ void indme(
delete[] dref;
}
-void instr(ScattererConfiguration *sconf, ParticleDescriptor *c1) {
+void instr(double **rcf, ParticleDescriptor *c1) {
const int ylm_size = (c1->litpos > c1->lmtpos) ? c1->litpos : c1->lmtpos;
dcomplex *ylm = new dcomplex[ylm_size]();
double rx, ry, rz, rr, crth, srth, crph, srph;
@@ -503,7 +503,7 @@ void instr(ScattererConfiguration *sconf, ParticleDescriptor *c1) {
if (c1->iog[i18] >= i) {
int nsh = c1->nshl[i18];
for (int j = 0; j < nsh; j++)
- c1->rc[i18][j] = sconf->get_rcf(i18, j) * c1->ros[i18];
+ c1->rc[i18][j] = rcf[i18][j] * c1->ros[i18];
}
} // i18 loop
int i = 0;
diff --git a/src/libnptm/outputs.cpp b/src/libnptm/outputs.cpp
index 48791dc9a947ef24f9183d974b5855d7ebe24602..8c4ecf7df21e3ed3d04320eeac8ec58666d2c510 100644
--- a/src/libnptm/outputs.cpp
+++ b/src/libnptm/outputs.cpp
@@ -4281,22 +4281,22 @@ int InclusionOutputInfo::write_legacy(const std::string &output) {
if (isam >= 0) {
fprintf(
p_outfile, " UNI=(%12.5lE,%12.5lE,%12.5lE)\n",
- vec_dir_un[done_dirs],
- vec_dir_un[done_dirs + 1],
- vec_dir_un[done_dirs + 2]
+ vec_dir_un[3 * done_dirs],
+ vec_dir_un[3 * done_dirs + 1],
+ vec_dir_un[3 * done_dirs + 2]
);
fprintf(
p_outfile, " UNS=(%12.5lE,%12.5lE,%12.5lE)\n",
- vec_dir_uns[done_dirs],
- vec_dir_uns[done_dirs + 1],
- vec_dir_uns[done_dirs + 2]
+ vec_dir_uns[3 * done_dirs],
+ vec_dir_uns[3 * done_dirs + 1],
+ vec_dir_uns[3 * done_dirs + 2]
);
} else { // label 214
fprintf(
p_outfile, " UN=(%12.5lE,%12.5lE,%12.5lE)\n\n",
- vec_dir_un[done_dirs],
- vec_dir_un[done_dirs + 1],
- vec_dir_un[done_dirs + 2]
+ vec_dir_un[3 * done_dirs],
+ vec_dir_un[3 * done_dirs + 1],
+ vec_dir_un[3 * done_dirs + 2]
);
}
fprintf(p_outfile, " SINGLE SCATTERER\n");
@@ -5644,7 +5644,7 @@ int SphereOutputInfo::write_legacy(const std::string &file_name) {
int done_dirs = 0;
fprintf(p_outfile, "========== JXI =%3d ====================\n", jxi + 1);
if (idfc >= 0) {
- fprintf(p_outfile, " VK=%15.7lE, XI=%15.7lE\n", vec_xi[jxi], vec_vk[jxi]);
+ fprintf(p_outfile, " VK=%15.7lE, XI=%15.7lE\n", vec_vk[jxi], vec_xi[jxi]);
} else { // IDFC < 0
fprintf(p_outfile, " XI=%15.7lE\n", vec_xi[jxi]);
}
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index 288c259f4aee6b48efe3b6980af1e87f01e6dbd8..7d26091f21762d4f1edf166e1a65cf6dd7c92163 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -201,25 +201,25 @@ void sphere(const string& config_file, const string& data_file, const string& ou
}
// Do the first wavelength iteration
- int jxi488 = 1;
+ // int jxi488 = 1;
// Use pragmas to put OMP parallelism to second level.
int jer = 0;
-#pragma omp parallel
- {
-#pragma omp single
- {
- jer = sphere_jxi488_cycle(jxi488 - 1, sconf, gconf, p_sa, sid, p_output, output_path, vtppoanp);
- } // OMP single
- } // OMP parallel
- if (jer != 0) { // First iteration failed. Halt the calculation.
- delete p_output;
- delete p_sa;
- delete sid;
- delete logger;
- delete sconf;
- delete gconf;
- return;
- }
+// #pragma omp parallel
+// {
+// #pragma omp single
+// {
+// jer = sphere_jxi488_cycle(jxi488 - 1, sconf, gconf, p_sa, sid, p_output, output_path, vtppoanp);
+// } // OMP single
+// } // OMP parallel
+// if (jer != 0) { // First iteration failed. Halt the calculation.
+// delete p_output;
+// delete p_sa;
+// delete sid;
+// delete logger;
+// delete sconf;
+// delete gconf;
+// return;
+// }
//==================================================
// do the first outputs here, so that I open here the new files, afterwards I only append
@@ -310,7 +310,7 @@ void sphere(const string& config_file, const string& data_file, const string& ou
}
#pragma omp barrier
// ok, now I can actually start the parallel calculations
- for (int ixi488=2; ixi488<=sid_2->number_of_scales; ixi488 +=myMPIstride) {
+ for (int ixi488=1; ixi488<=sid_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;
@@ -463,7 +463,7 @@ void sphere(const string& config_file, const string& data_file, const string& ou
// 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<=sid_2->number_of_scales; ixi488 +=myMPIstride) {
+ for (int ixi488=1; ixi488<=sid_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;
@@ -542,11 +542,15 @@ 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;
+ // this is now part of sid, so don't mess with it here, just copy it by reference
+ bool &is_first_scale = sid->is_first_scale;
+ //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,6 +580,30 @@ int sphere_jxi488_cycle(
oi->vec_vk[jxindex] = vk;
oi->vec_xi[jxindex] = xi;
}
+ // 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++) {
@@ -595,7 +623,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);
}
@@ -718,7 +746,7 @@ int sphere_jxi488_cycle(
double ph = sa->ph;
for (int jph484 = 0; jph484 < sa->nph; jph484++) {
int jph = jph484 + 1;
- bool goto182 = (sa->nk == 1) && (jxi > 1);
+ bool goto182 = (sa->nk == 1) && (!is_first_scale);
if (!goto182) {
upvmp(th, ph, 0, sid->cost, sid->sint, sid->cosp, sid->sinp, sid->u, sid->upmp, sid->unmp);
}
@@ -753,14 +781,14 @@ int sphere_jxi488_cycle(
phs = ph + phsph;
if (phs >= 360.0) phs -= 360.0;
}
- bool goto190 = (nks == 1) && ((jxi > 1) || (jth > 1) || (jph > 1));
+ bool goto190 = (nks == 1) && ((!is_first_scale) || (jth > 1) || (jph > 1));
if (!goto190) {
upvmp(ths, phs, icspnv, sid->costs, sid->sints, sid->cosps, sid->sinps, sid->us, sid->upsmp, sid->unsmp);
if (gconf->isam >= 0) {
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,
@@ -863,6 +891,8 @@ int sphere_jxi488_cycle(
} // jph484 loop on elevation
th += sa->thstp;
} // jth486 loop on azimuth
+ // 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;
oi->vec_jxi[jxindex] = jxi;
logger->log("INFO: finished scale iteration " + to_string(jxi) + " of " + to_string(nxi) + ".\n");
delete logger;
@@ -951,6 +981,7 @@ SphereIterationData::SphereIterationData(
zpv[zi][zj][1] = (vec_zpv + vec_index + 2);
}
}
+ is_first_scale = 1;
}
SphereIterationData::SphereIterationData(const SphereIterationData &rhs) {
@@ -1056,6 +1087,7 @@ SphereIterationData::SphereIterationData(const SphereIterationData &rhs) {
zpv[zi][zj][1] = (vec_zpv + vec_index + 2);
}
}
+ is_first_scale = rhs.is_first_scale;
}
SphereIterationData::~SphereIterationData() {
@@ -1181,6 +1213,9 @@ SphereIterationData::SphereIterationData(const mixMPI *mpidata, const int device
// Collect vectors whose size depend on LM
vec_zpv = new double[12 * _lm];
MPI_Bcast(vec_zpv, 12 * _lm, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+
+ MPI_Bcast(&is_first_scale, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);
+
}
int SphereIterationData::mpibcast(const mixMPI *mpidata) {
@@ -1243,8 +1278,40 @@ int SphereIterationData::mpibcast(const mixMPI *mpidata) {
// Broadcast vectors whose size depend on LM
MPI_Bcast(vec_zpv, 12 * _lm, MPI_DOUBLE, 0, MPI_COMM_WORLD);
-
+
+ MPI_Bcast(&is_first_scale, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);
+
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 <<<