From 71e7f35be62b7845cfa3b6834e18ba6138d22c7d Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Tue, 18 Mar 2025 15:54:08 +0100
Subject: [PATCH] Adapt dc0 vector in ParticlDescriptor for multi-layer
aggregate case
---
src/include/Commons.h | 8 ++++----
src/libnptm/Commons.cpp | 30 ++++++++++++++++++------------
2 files changed, 22 insertions(+), 16 deletions(-)
diff --git a/src/include/Commons.h b/src/include/Commons.h
index cce0a653..31efcc39 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -80,12 +80,12 @@ protected:
int _nsph;
//! \brief Maximum internal field expansion order.
int _li;
+ //! \brief Maximum number of layers in known sphere types.
+ int _max_layers;
//! \brief Number of different sphere types.
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.
@@ -184,12 +184,12 @@ public:
const int &nsph = _nsph;
//! \brief Read-only view of maximum internal field expansion order.
const int &li = _li;
+ //! \brief Read-only view of the maximum number of layers in types.
+ const int &max_layers = _max_layers;
//! \brief Read-only view of number of different sphere types.
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.
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index cfdec141..4a29057c 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -84,7 +84,12 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo
gcs = 0.0;
_num_configurations = sconf->configurations;
_num_layers = (sconf->use_external_sphere) ? 1 : 0;
- for (int nli = 0; nli < num_configurations; nli++) _num_layers += sconf->get_nshl(nli);
+ _max_layers = 1;
+ for (int nli = 0; nli < num_configurations; nli++) {
+ int nl = sconf->get_nshl(nli);
+ _num_layers += nl;
+ if (nl > _max_layers) _max_layers = nl;
+ }
vec_rmi = new dcomplex[_li * _nsph]();
rmi = new dcomplex*[_li];
@@ -131,12 +136,12 @@ 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;
+ // _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]();
+ dc0 = new dcomplex[_max_layers + 1]();
vsz = new double[_nsph]();
// >>> NEEDED BY SPHERE AND CLUSTER <<<
@@ -219,6 +224,7 @@ ParticleDescriptor::ParticleDescriptor(const ParticleDescriptor &rhs) {
_ndit = rhs._ndit;
_ndm = rhs._ndm;
gcs = rhs.gcs;
+ _max_layers = rhs._max_layers;
_num_configurations = rhs._num_configurations;
_num_layers = rhs._num_layers;
@@ -271,7 +277,7 @@ ParticleDescriptor::ParticleDescriptor(const ParticleDescriptor &rhs) {
_npnt = rhs._npnt;
_npntts = rhs._npntts;
_nhspo = rhs._nhspo;
- _nl = rhs._nl;
+ _max_layers = rhs._max_layers;
ris = new dcomplex[_nhspo]();
dlri = new dcomplex[_nhspo]();
for (int ri = 0; ri < _nhspo; ri++) {
@@ -284,8 +290,8 @@ ParticleDescriptor::ParticleDescriptor(const ParticleDescriptor &rhs) {
vkt[vi] = rhs.vkt[vi];
vsz[vi] = rhs.vsz[vi];
}
- dc0 = new dcomplex[_nl]();
- for (int di = 0; di < _nl; di++) {
+ dc0 = new dcomplex[_max_layers]();
+ for (int di = 0; di < _max_layers; di++) {
dc0[di] = rhs.dc0[di];
}
// >>> NEEDED BY SPHERE AND CLUSTER <<<
@@ -416,7 +422,7 @@ ParticleDescriptor::ParticleDescriptor(const mixMPI *mpidata) {
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(&_max_layers, 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];
@@ -425,8 +431,8 @@ ParticleDescriptor::ParticleDescriptor(const mixMPI *mpidata) {
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);
+ dc0 = new dcomplex[_max_layers];
+ MPI_Bcast(dc0, _max_layers, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
}
void ParticleDescriptor::mpibcast(const mixMPI *mpidata) {
@@ -466,12 +472,12 @@ void ParticleDescriptor::mpibcast(const mixMPI *mpidata) {
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(&_max_layers, 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);
+ MPI_Bcast(dc0, _max_layers, 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);
--
GitLab