Implement CLUSTER dynamic handling of internal order

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>
@@ -240,7 +241,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);
@@ -680,6 +681,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
 {
   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 +691,24 @@ 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;
+  // int li = gconf->li;
+  // int le = gconf->le;
+  // int lm = 0;
+  // if (le > lm) lm = le;
+  // if (li > lm) lm = li;
+  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;
+  // np_int ndit = 2 * nsph * cid->c1->nlim;
   int isq, ibf;
   int last_configuration;
   int num_configs = sconf->configurations;
-  int ndirs = sa->nkks;
+  const int ndirs = sa->nkks;
   int oindex = -1;
   int jindex = jxi488 - output->first_xi + 1;
 
@@ -716,7 +718,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 +731,29 @@ 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 double alamb = 2.0 * pi / cid->vk;
+  double size_par = 2.0 * pi * sqrt(exdc) * sconf->get_radius(0) / alamb;
+  int recommended_li = 4 + (int)ceil(size_par + 4.05 * pow(size_par, 1.0 / 3.0));
+  if (recommended_li != cid->c1->li) {
+    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 {
+      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);
+      cid->update_orders(sconf->_rcf, recommended_li, recommended_li);
+    }
+    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;
@@ -837,6 +862,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";
@@ -2062,4 +2088,47 @@ ClusterIterationData::~ClusterIterationData() {
   delete[] cmullr;
   delete[] cmul;
 }
+
+int ClusterIterationData::update_orders(double **rcf, int inner_order, int outer_order) {
+  int result = 0;
+  int old_li = c1->li;
+  int old_le = c1->le;
+  int old_lm = c1->lm;
+  np_int old_ndit = 2 * c1->nsph * c1->nlim;
+  if (inner_order != old_li || outer_order != old_le) {
+    ((ParticleDescriptorCluster *)c1)->update_orders(inner_order, outer_order);
+    const int ndi = c1->nsph * c1->nlim;
+    const np_int ndit = 2 * ndi;
+    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 <<<