diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp index f3e9141dd5567e96a0155c03099870a1f6d7fcd8..7fa92b0eb9241cd92550cf310e3c357176f6981a 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> @@ -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 <<<