Skip to content
Snippets Groups Projects
Commit 7148fba0 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Use InclusionOutputInfo in np_inclusion

parent 34cbccdb
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/inclusion/inclusion.cpp
+ 250
− 279
View file @ 7148fba0
...@@ -91,6 +91,10 @@ ...@@ -91,6 +91,10 @@
#include "../include/utils.h" #include "../include/utils.h"
#endif #endif
#ifndef INCLUDE_OUTPUTS_H_
#include "../include/outputs.h"
#endif
using namespace std; using namespace std;
// >>> InclusionIterationData header <<< // // >>> InclusionIterationData header <<< //
...@@ -816,11 +820,11 @@ InclusionIterationData::~InclusionIterationData() { ...@@ -816,11 +820,11 @@ InclusionIterationData::~InclusionIterationData() {
* \param gconf: `GeometryConfiguration *` Pointer to a `GeometryConfiguration` object. * \param gconf: `GeometryConfiguration *` Pointer to a `GeometryConfiguration` object.
* \param sa: `ScatteringAngles *` Pointer to a `ScatteringAngles` object. * \param sa: `ScatteringAngles *` Pointer to a `ScatteringAngles` object.
* \param cid: `InclusionIterationData *` Pointer to an `InclusionIterationData` object. * \param cid: `InclusionIterationData *` Pointer to an `InclusionIterationData` object.
* \param output: `VirtualAsciiFile *` Pointer to a `VirtualAsciiFile` object. * \param oi: `InclusionOutputInfo *` Pointer to an `InclusionOutputInfo` object.
* \param output_path: `const string &` Path to the output directory. * \param output_path: `const string &` Path to the output directory.
* \param vtppoanp: `VirtualBinaryFile *` Pointer to a `VirtualBinaryFile` object. * \param vtppoanp: `VirtualBinaryFile *` Pointer to a `VirtualBinaryFile` object.
*/ */
int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, InclusionIterationData *cid, VirtualAsciiFile *output, const string& output_path, VirtualBinaryFile *vtppoanp); int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, InclusionIterationData *cid, InclusionOutputInfo *output, const string& output_path, VirtualBinaryFile *vtppoanp);
/*! \brief C++ implementation of INCLU /*! \brief C++ implementation of INCLU
* *
...@@ -929,81 +933,17 @@ void inclusion(const string& config_file, const string& data_file, const string& ...@@ -929,81 +933,17 @@ void inclusion(const string& config_file, const string& data_file, const string&
ScatteringAngles *p_scattering_angles = new ScatteringAngles(gconf); ScatteringAngles *p_scattering_angles = new ScatteringAngles(gconf);
double wp = sconf->wp; double wp = sconf->wp;
// Open empty virtual ascii file for output // Open empty virtual ascii file for output
VirtualAsciiFile *p_output = new VirtualAsciiFile(); InclusionOutputInfo *p_output = new InclusionOutputInfo(sconf, gconf, mpidata);
char virtual_line[256];
InclusionIterationData *cid = new InclusionIterationData(gconf, sconf, mpidata, device_count); InclusionIterationData *cid = new InclusionIterationData(gconf, sconf, mpidata, device_count);
const np_int ndi = cid->c1->nsph * cid->c1->nlim; const np_int ndi = cid->c1->nsph * cid->c1->nlim;
const np_int ndit = 2 * ndi; 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"); 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"); 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");
//==========================
// Write a block of info to the ascii output file
//==========================
sprintf(virtual_line, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
p_output->append_line(virtual_line);
#ifdef USE_ILP64
sprintf(virtual_line, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",
nsph, cid->c1->li, cid->c1->le, gconf->mxndm, gconf->in_pol, gconf->npnt,
gconf->npntts, gconf->iavm, gconf->iavm
);
#else
sprintf(virtual_line, " %5d%5d%5d%5d%5d%5d%5d%5d%5d\n",
nsph, cid->c1->li, cid->c1->le, gconf->mxndm, gconf->in_pol, gconf->npnt,
gconf->npntts, gconf->iavm, gconf->iavm
);
#endif // USE_ILP64
p_output->append_line(virtual_line);
sprintf(virtual_line, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
p_output->append_line(virtual_line);
for (int ri = 0; ri < nsph; ri++) {
sprintf(virtual_line, "%17.8lE%17.8lE%17.8lE\n",
gconf->get_sph_x(ri), gconf->get_sph_y(ri), gconf->get_sph_z(ri)
);
p_output->append_line(virtual_line);
}
sprintf(virtual_line, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
p_output->append_line(virtual_line);
sprintf(
virtual_line, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
p_scattering_angles->th, p_scattering_angles->thstp,
p_scattering_angles->thlst, p_scattering_angles->ths,
p_scattering_angles->thsstp, p_scattering_angles->thslst
);
p_output->append_line(virtual_line);
sprintf(virtual_line, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
p_output->append_line(virtual_line);
sprintf(
virtual_line, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
p_scattering_angles->ph, p_scattering_angles->phstp,
p_scattering_angles->phlst, p_scattering_angles->phs,
p_scattering_angles->phsstp, p_scattering_angles->phslst
);
p_output->append_line(virtual_line);
sprintf(virtual_line, " READ(IR,*)JWTM\n");
p_output->append_line(virtual_line);
sprintf(virtual_line, " %5d\n", gconf->jwtm);
p_output->append_line(virtual_line);
sprintf(virtual_line, " READ(ITIN)NSPHT\n");
p_output->append_line(virtual_line);
sprintf(virtual_line, " READ(ITIN)(IOG(I),I=1,NSPH)\n");
p_output->append_line(virtual_line);
sprintf(virtual_line, " READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");
p_output->append_line(virtual_line);
sprintf(virtual_line, " READ(ITIN)(XIV(I),I=1,NXI)\n");
p_output->append_line(virtual_line);
sprintf(virtual_line, " READ(ITIN)NSHL(I),ROS(I)\n");
p_output->append_line(virtual_line);
sprintf(virtual_line, " READ(ITIN)(RCF(I,NS),NS=1,NSH)\n");
p_output->append_line(virtual_line);
sprintf(virtual_line, " \n");
p_output->append_line(virtual_line);
instr(sconf, cid->c1); instr(sconf, cid->c1);
thdps(cid->c1->lm, cid->zpv); thdps(cid->c1->lm, cid->zpv);
double exdc = sconf->exdc; double exdc = sconf->exdc;
double exri = sqrt(exdc); double exri = sqrt(exdc);
sprintf(virtual_line, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
p_output->append_line(virtual_line);
// Create an empty bynary file // Create an empty bynary file
VirtualBinaryFile *vtppoanp = new VirtualBinaryFile(); VirtualBinaryFile *vtppoanp = new VirtualBinaryFile();
...@@ -1037,10 +977,7 @@ void inclusion(const string& config_file, const string& data_file, const string& ...@@ -1037,10 +977,7 @@ void inclusion(const string& config_file, const string& data_file, const string&
vtppoanp->append_line(VirtualBinaryLine(nphs)); vtppoanp->append_line(VirtualBinaryLine(nphs));
if (sconf->idfc < 0) { if (sconf->idfc < 0) {
cid->vk = cid->xip * cid->wn; cid->vk = cid->xip * cid->wn;
sprintf(virtual_line, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", cid->vk); p_output->vec_vk[0] = cid->vk;
p_output->append_line(virtual_line);
sprintf(virtual_line, " \n");
p_output->append_line(virtual_line);
} }
// do the first iteration on jxi488 separately, since it seems to be different from the others // do the first iteration on jxi488 separately, since it seems to be different from the others
...@@ -1092,8 +1029,6 @@ void inclusion(const string& config_file, const string& data_file, const string& ...@@ -1092,8 +1029,6 @@ void inclusion(const string& config_file, const string& data_file, const string&
//================================================== //==================================================
// do the first outputs here, so that I open here the new files, afterwards I only append // do the first outputs here, so that I open here the new files, afterwards I only append
//================================================== //==================================================
p_output->write_to_disk(output_path + "/c_OINCLU");
delete p_output;
vtppoanp->write_to_disk(output_path + "/c_TPPOAN"); vtppoanp->write_to_disk(output_path + "/c_TPPOAN");
delete vtppoanp; delete vtppoanp;
...@@ -1113,7 +1048,7 @@ void inclusion(const string& config_file, const string& data_file, const string& ...@@ -1113,7 +1048,7 @@ void inclusion(const string& config_file, const string& data_file, const string&
int myMPIstride = ompnumthreads; int myMPIstride = ompnumthreads;
int myMPIblock = ompnumthreads; int myMPIblock = ompnumthreads;
// Define here shared arrays of virtual ascii and binary files, so that thread 0 will be able to access them all later // Define here shared arrays of virtual ascii and binary files, so that thread 0 will be able to access them all later
VirtualAsciiFile **p_outarray = NULL; InclusionOutputInfo **p_outarray = NULL;
VirtualBinaryFile **vtppoanarray = NULL; VirtualBinaryFile **vtppoanarray = NULL;
#ifdef USE_NVTX #ifdef USE_NVTX
...@@ -1136,7 +1071,7 @@ void inclusion(const string& config_file, const string& data_file, const string& ...@@ -1136,7 +1071,7 @@ void inclusion(const string& config_file, const string& data_file, const string&
if (myompthread == 0) { if (myompthread == 0) {
// Initialise some shared variables only on thread 0 // Initialise some shared variables only on thread 0
p_outarray = new VirtualAsciiFile*[ompnumthreads]; p_outarray = new InclusionOutputInfo*[ompnumthreads];
vtppoanarray = new VirtualBinaryFile*[ompnumthreads]; vtppoanarray = new VirtualBinaryFile*[ompnumthreads];
myMPIblock = ompnumthreads; myMPIblock = ompnumthreads;
myMPIstride = myMPIblock; myMPIstride = myMPIblock;
...@@ -1165,11 +1100,14 @@ void inclusion(const string& config_file, const string& data_file, const string& ...@@ -1165,11 +1100,14 @@ void inclusion(const string& config_file, const string& data_file, const string&
// To test parallelism, I will now start feeding this function with "clean" copies of the parameters, so that they will not be changed by previous iterations, and each one will behave as the first one. Define all (empty) variables here, so they have the correct scope, then they get different definitions depending on thread number // To test parallelism, I will now start feeding this function with "clean" copies of the parameters, so that they will not be changed by previous iterations, and each one will behave as the first one. Define all (empty) variables here, so they have the correct scope, then they get different definitions depending on thread number
InclusionIterationData *cid_2 = NULL; InclusionIterationData *cid_2 = NULL;
VirtualAsciiFile *p_output_2 = NULL; InclusionOutputInfo *p_output_2 = NULL;
VirtualBinaryFile *vtppoanp_2 = NULL; VirtualBinaryFile *vtppoanp_2 = NULL;
// for threads other than the 0, create distinct copies of all relevant data, while for thread 0 just define new references / pointers to the original ones // for threads other than the 0, create distinct copies of all relevant data, while for thread 0 just define new references / pointers to the original ones
if (myompthread == 0) { if (myompthread == 0) {
cid_2 = cid; cid_2 = cid;
// OMP thread 0 of MPI process 0 holds the pointer to the full output structure
p_output_2 = p_output;
p_outarray[0] = p_output_2;
} else { } else {
// this is not thread 0, so do create fresh copies of all local variables // this is not thread 0, so do create fresh copies of all local variables
cid_2 = new InclusionIterationData(*cid); cid_2 = new InclusionIterationData(*cid);
...@@ -1185,16 +1123,24 @@ void inclusion(const string& config_file, const string& data_file, const string& ...@@ -1185,16 +1123,24 @@ void inclusion(const string& config_file, const string& data_file, const string&
#pragma omp barrier #pragma omp barrier
int myjxi488 = ixi488+myompthread; int myjxi488 = ixi488+myompthread;
// each thread opens new virtual files and stores their pointers in the shared array // each thread opens new virtual files and stores their pointers in the shared array
p_output_2 = new VirtualAsciiFile();
vtppoanp_2 = new VirtualBinaryFile(); vtppoanp_2 = new VirtualBinaryFile();
// each thread puts a copy of the pointers to its virtual files in the shared arrays // each thread puts a copy of the pointers to its virtual files in the shared arrays
p_outarray[myompthread] = p_output_2;
vtppoanarray[myompthread] = vtppoanp_2; vtppoanarray[myompthread] = vtppoanp_2;
#pragma omp barrier #pragma omp barrier
// each MPI process handles a number of contiguous scales corresponding to its number of OMP threads at this omp level of parallelism // each MPI process handles a number of contiguous scales corresponding to its number of OMP threads at this omp level of parallelism
if (myjxi488 <= cid_2->number_of_scales) { if (myjxi488 <= cid_2->number_of_scales) {
if (myompthread > 0) {
// UPDATE: non-0 threads need to allocate memory for one scale at a time.
p_output_2 = new InclusionOutputInfo(sconf, gconf, mpidata, myjxi488, 1);
p_outarray[myompthread] = p_output_2;
}
int jer = inclusion_jxi488_cycle(myjxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2); int jer = inclusion_jxi488_cycle(myjxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2);
} else {
if (myompthread > 0) {
// If there is no input for this thread, set output pointer to NULL.
p_outarray[myompthread] = NULL;
}
} }
#pragma omp barrier #pragma omp barrier
...@@ -1205,8 +1151,11 @@ void inclusion(const string& config_file, const string& data_file, const string& ...@@ -1205,8 +1151,11 @@ void inclusion(const string& config_file, const string& data_file, const string&
// threads different from 0 append their virtual files to the one of thread 0, and delete them // threads different from 0 append their virtual files to the one of thread 0, and delete them
if (myompthread == 0) { if (myompthread == 0) {
for (int ti=1; ti<ompnumthreads; ti++) { for (int ti=1; ti<ompnumthreads; ti++) {
p_outarray[0]->append(*(p_outarray[ti])); if (p_outarray[ti] != NULL) {
p_outarray[0]->insert(*(p_outarray[ti]));
delete p_outarray[ti]; delete p_outarray[ti];
p_outarray[ti] = NULL;
}
vtppoanarray[0]->append(*(vtppoanarray[ti])); vtppoanarray[0]->append(*(vtppoanarray[ti]));
delete vtppoanarray[ti]; delete vtppoanarray[ti];
} }
...@@ -1217,8 +1166,8 @@ void inclusion(const string& config_file, const string& data_file, const string& ...@@ -1217,8 +1166,8 @@ void inclusion(const string& config_file, const string& data_file, const string&
//============================================== //==============================================
if (myompthread == 0) { if (myompthread == 0) {
// thread 0 writes its virtual files, now including contributions from all threads, to disk, and deletes them // 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_OINCLU"); // p_outarray[0]->append_to_disk(output_path + "/c_OINCLU");
delete p_outarray[0]; // delete p_outarray[0];
vtppoanarray[0]->append_to_disk(output_path + "/c_TPPOAN"); vtppoanarray[0]->append_to_disk(output_path + "/c_TPPOAN");
delete vtppoanarray[0]; delete vtppoanarray[0];
...@@ -1226,11 +1175,14 @@ void inclusion(const string& config_file, const string& data_file, const string& ...@@ -1226,11 +1175,14 @@ void inclusion(const string& config_file, const string& data_file, const string&
if (mpidata->mpirunning) { if (mpidata->mpirunning) {
// only go through this if MPI has been actually used // only go through this if MPI has been actually used
for (int rr=1; rr<mpidata->nprocs; rr++) { for (int rr=1; rr<mpidata->nprocs; rr++) {
// get the data from process rr by receiving it in total memory structure
p_outarray[0]->mpireceive(mpidata, rr);
// get the data from process rr, creating a new virtual ascii file // get the data from process rr, creating a new virtual ascii file
VirtualAsciiFile *p_output = new VirtualAsciiFile(mpidata, rr); // VirtualAsciiFile *p_output = new VirtualAsciiFile(mpidata, rr);
// append to disk and delete virtual ascii file // append to disk and delete virtual ascii file
p_output->append_to_disk(output_path + "/c_OINCLU"); // p_output->append_to_disk(output_path + "/c_OINCLU");
delete p_output; // delete p_output;
// get the data from process rr, creating a new virtual binary file // get the data from process rr, creating a new virtual binary file
VirtualBinaryFile *vtppoanp = new VirtualBinaryFile(mpidata, rr); VirtualBinaryFile *vtppoanp = new VirtualBinaryFile(mpidata, rr);
// append to disk and delete virtual binary file // append to disk and delete virtual binary file
...@@ -1264,7 +1216,10 @@ void inclusion(const string& config_file, const string& data_file, const string& ...@@ -1264,7 +1216,10 @@ void inclusion(const string& config_file, const string& data_file, const string&
delete cid_2; delete cid_2;
} }
delete p_scattering_angles; delete p_scattering_angles;
} p_output->write(output_path + "/c_OINCLU.hd5", "HDF5");
p_output->write(output_path + "/c_OINCLU", "LEGACY");
delete p_output;
} // closes s_nsph == nsph check
else { // Sphere number inconsistency error. else { // Sphere number inconsistency error.
throw UnrecognizedConfigurationException( throw UnrecognizedConfigurationException(
...@@ -1302,7 +1257,7 @@ void inclusion(const string& config_file, const string& data_file, const string& ...@@ -1302,7 +1257,7 @@ void inclusion(const string& config_file, const string& data_file, const string&
// Create this variable and initialise it with a default here, so that it is defined anyway, with or without OpenMP support enabled // Create this variable and initialise it with a default here, so that it is defined anyway, with or without OpenMP support enabled
int ompnumthreads = 1; int ompnumthreads = 1;
VirtualAsciiFile **p_outarray = NULL; InclusionOutputInfo **p_outarray = NULL;
VirtualBinaryFile **vtppoanarray = NULL; VirtualBinaryFile **vtppoanarray = NULL;
int myjxi488startoffset; int myjxi488startoffset;
int myMPIstride = ompnumthreads; int myMPIstride = ompnumthreads;
...@@ -1325,13 +1280,13 @@ void inclusion(const string& config_file, const string& data_file, const string& ...@@ -1325,13 +1280,13 @@ void inclusion(const string& config_file, const string& data_file, const string&
// receive myMPIstride sent by MPI process 0 to all processes // receive myMPIstride sent by MPI process 0 to all processes
MPI_Bcast(&myMPIstride, 1, MPI_INT, 0, MPI_COMM_WORLD); MPI_Bcast(&myMPIstride, 1, MPI_INT, 0, MPI_COMM_WORLD);
// allocate virtual files for each thread // allocate virtual files for each thread
p_outarray = new VirtualAsciiFile*[ompnumthreads]; p_outarray = new InclusionOutputInfo*[ompnumthreads];
vtppoanarray = new VirtualBinaryFile*[ompnumthreads]; vtppoanarray = new VirtualBinaryFile*[ompnumthreads];
} }
#pragma omp barrier #pragma omp barrier
// To test parallelism, I will now start feeding this function with "clean" copies of the parameters, so that they will not be changed by previous iterations, and each one will behave as the first one. Define all (empty) variables here, so they have the correct scope, then they get different definitions depending on thread number // To test parallelism, I will now start feeding this function with "clean" copies of the parameters, so that they will not be changed by previous iterations, and each one will behave as the first one. Define all (empty) variables here, so they have the correct scope, then they get different definitions depending on thread number
InclusionIterationData *cid_2 = NULL; InclusionIterationData *cid_2 = NULL;
VirtualAsciiFile *p_output_2 = NULL; InclusionOutputInfo *p_output_2 = NULL;
VirtualBinaryFile *vtppoanp_2 = NULL; VirtualBinaryFile *vtppoanp_2 = NULL;
// PLACEHOLDER // PLACEHOLDER
// for threads other than the 0, create distinct copies of all relevant data, while for thread 0 just define new references / pointers to the original ones // for threads other than the 0, create distinct copies of all relevant data, while for thread 0 just define new references / pointers to the original ones
...@@ -1349,25 +1304,41 @@ void inclusion(const string& config_file, const string& data_file, const string& ...@@ -1349,25 +1304,41 @@ void inclusion(const string& config_file, const string& data_file, const string&
#pragma omp barrier #pragma omp barrier
int myjxi488 = ixi488 + myjxi488startoffset + myompthread; int myjxi488 = ixi488 + myjxi488startoffset + myompthread;
// each thread opens new virtual files and stores their pointers in the shared array // each thread opens new virtual files and stores their pointers in the shared array
p_output_2 = new VirtualAsciiFile();
vtppoanp_2 = new VirtualBinaryFile(); vtppoanp_2 = new VirtualBinaryFile();
// each thread puts a copy of the pointers to its virtual files in the shared arrays // each thread puts a copy of the pointers to its virtual files in the shared arrays
p_outarray[myompthread] = p_output_2;
vtppoanarray[myompthread] = vtppoanp_2; vtppoanarray[myompthread] = vtppoanp_2;
#pragma omp barrier #pragma omp barrier
if (myompthread==0) logger->log("Syncing OpenMP threads and starting the loop on wavelengths\n"); if (myompthread==0) logger->log("Syncing OpenMP threads and starting the loop on wavelengths\n");
// ok, now I can actually start the parallel calculations // ok, now I can actually start the parallel calculations
// each MPI process handles a number of contiguous scales corresponding to its number of OMP threads at this omp level of parallelism // each MPI process handles a number of contiguous scales corresponding to its number of OMP threads at this omp level of parallelism
if (myjxi488 <= cid_2->number_of_scales) { if (myjxi488 <= cid_2->number_of_scales) {
if (myompthread > 0) {
// UPDATE: non-0 threads need to allocate memory for one scale at a time.
p_output_2 = new InclusionOutputInfo(sconf, gconf, mpidata, myjxi488, 1);
p_outarray[myompthread] = p_output_2;
} else {
// Thread 0 of non-zero MPI processes needs to allocate memory for the
// output of all threads.
p_output_2 = new InclusionOutputInfo(sconf, gconf, mpidata, myjxi488, ompnumthreads);
p_outarray[0] = p_output_2;
}
int jer = inclusion_jxi488_cycle(myjxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2); int jer = inclusion_jxi488_cycle(myjxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2);
} // close the OMP parallel for loop } else {
if (myompthread > 0) {
// If there is no input for this thread, set the output pointer to NULL.
p_outarray[myompthread] = NULL;
}
}
#pragma omp barrier #pragma omp barrier
// threads different from 0 append their virtual files to the one of thread 0, and delete them // threads different from 0 append their virtual files to the one of thread 0, and delete them
if (myompthread == 0) { if (myompthread == 0) {
for (int ti=1; ti<ompnumthreads; ti++) { for (int ti=1; ti<ompnumthreads; ti++) {
p_outarray[0]->append(*(p_outarray[ti])); if (p_outarray[ti] != NULL) {
p_outarray[0]->insert(*(p_outarray[ti]));
delete p_outarray[ti]; delete p_outarray[ti];
p_outarray[ti] = NULL;
}
vtppoanarray[0]->append(*(vtppoanarray[ti])); vtppoanarray[0]->append(*(vtppoanarray[ti]));
delete vtppoanarray[ti]; delete vtppoanarray[ti];
} }
...@@ -1407,7 +1378,8 @@ void inclusion(const string& config_file, const string& data_file, const string& ...@@ -1407,7 +1378,8 @@ void inclusion(const string& config_file, const string& data_file, const string&
#endif #endif
} }
int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, InclusionIterationData *cid, VirtualAsciiFile *output, const string& output_path, VirtualBinaryFile *vtppoanp) { 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; int nxi = sconf->number_of_scales;
char virtual_line[256]; char virtual_line[256];
string message = "INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n"; string message = "INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n";
...@@ -1434,12 +1406,14 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo ...@@ -1434,12 +1406,14 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
int last_configuration; int last_configuration;
dcomplex ent, entn; dcomplex ent, entn;
double enti; double enti;
int num_configs = sconf->configurations;
int ndirs = sa->nkks;
int oindex = -1;
int jindex = jxi488 - output->first_xi + 1;
#ifdef USE_NVTX #ifdef USE_NVTX
nvtxRangePush("Prepare matrix calculation"); nvtxRangePush("Prepare matrix calculation");
#endif #endif
sprintf(virtual_line, "========== JXI =%3d ====================\n", jxi488);
output->append_line(virtual_line);
double xi = sconf->get_scale(jxi488 - 1); double xi = sconf->get_scale(jxi488 - 1);
double exdc = sconf->exdc; double exdc = sconf->exdc;
double exri = sqrt(exdc); double exri = sqrt(exdc);
...@@ -1448,14 +1422,14 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo ...@@ -1448,14 +1422,14 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
if (idfc >= 0) { if (idfc >= 0) {
cid->vk = xi * cid->wn; cid->vk = xi * cid->wn;
vkarg = cid->vk; vkarg = cid->vk;
sprintf(virtual_line, " VK=%15.7lE, XI=%15.7lE\n", cid->vk, xi); output->vec_vk[jindex - 1] = cid->vk;
output->append_line(virtual_line); output->vec_xi[jindex - 1] = xi;
// goes to 120 // goes to 120
} else { // label 119 } else { // label 119
vkarg = xi * cid->vk; vkarg = xi * cid->vk;
cid->sqsfi = 1.0 / (xi * xi); cid->sqsfi = 1.0 / (xi * xi);
sprintf(virtual_line, " XI=%15.7lE\n", xi); output->vec_vk[jindex - 1] = cid->vk;
output->append_line(virtual_line); output->vec_xi[jindex - 1] = xi;
} }
// label 120 // label 120
double sze = vkarg * cid->extr; double sze = vkarg * cid->extr;
...@@ -1494,8 +1468,8 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo ...@@ -1494,8 +1468,8 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
} }
indme(i133, npnt, npntts, vkarg, ent, enti, entn, jer, lcalc, cid->arg, cid->c1); indme(i133, npnt, npntts, vkarg, ent, enti, entn, jer, lcalc, cid->arg, cid->c1);
if (jer != 0) { if (jer != 0) {
sprintf(virtual_line, " STOP IN INDME\n"); output->vec_ier[jindex - 1] = 1;
output->append_line(virtual_line); output->vec_jxi[jindex - 1] = -jxi488;
message = "ERROR: indme failed with error code " + to_string(jer) + ".\n"; message = "ERROR: indme failed with error code " + to_string(jer) + ".\n";
logger->log(message, LOG_ERRO); logger->log(message, LOG_ERRO);
delete logger; delete logger;
...@@ -1506,8 +1480,8 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo ...@@ -1506,8 +1480,8 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
} // i133 loop } // i133 loop
ospv(cid->c1, vkarg, sze, exri, entn, enti, jer, lcalc, cid->arg); ospv(cid->c1, vkarg, sze, exri, entn, enti, jer, lcalc, cid->arg);
if (jer != 0) { if (jer != 0) {
sprintf(virtual_line, " STOP IN OSPV\n"); output->vec_ier[jindex - 1] = 2;
output->append_line(virtual_line); output->vec_jxi[jindex - 1] = -jxi488;
message = "ERROR: ospv failed with error code " + to_string(jer) + ".\n"; message = "ERROR: ospv failed with error code " + to_string(jer) + ".\n";
logger->log(message, LOG_ERRO); logger->log(message, LOG_ERRO);
delete logger; delete logger;
...@@ -1610,19 +1584,25 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo ...@@ -1610,19 +1584,25 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
} }
} }
// label 156: continue from here // label 156: continue from here
last_configuration = 0;
for (int i168 = 1; i168 <= nsph; i168++) { for (int i168 = 1; i168 <= nsph; i168++) {
if (cid->c1->iog[i168 - 1] >= i168) { if (cid->c1->iog[i168 - 1] >= i168) {
int i = i168 - 1;
last_configuration++;
oindex = (jindex - 1) * (num_configs + 1) + last_configuration - 1;
if (cid->c1->nshl[i168 - 1] != 1) { if (cid->c1->nshl[i168 - 1] != 1) {
sprintf(virtual_line, " SPHERE N.%2d: SIZE=%15.7lE\n", i168, cid->c1->vsz[i168 - 1]); output->vec_sphere_ref_indices[oindex] = cc0;
output->append_line(virtual_line); output->vec_sphere_sizes[oindex] = cid->c1->vsz[i];
} else { } else {
sprintf(virtual_line, " SPHERE N.%2d: SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", i168, cid->c1->vsz[i168 - 1], real(cid->c1->vkt[i168 - 1]), imag(cid->c1->vkt[i168 - 1])); output->vec_sphere_ref_indices[oindex] = cid->c1->vkt[i];
output->append_line(virtual_line); output->vec_sphere_sizes[oindex] = cid->c1->vsz[i];
} }
} }
} // i168 loop } // i168 loop
sprintf(virtual_line, " EXT. SPHERE: SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", sze, real(entn), imag(entn)); last_configuration++;
output->append_line(virtual_line); oindex = (jindex - 1) * (num_configs + 1) + last_configuration - 1;
output->vec_sphere_sizes[oindex] = sze;
output->vec_sphere_ref_indices[oindex] = entn;
// label 160 // label 160
double cs0 = 0.25 * cid->vk * cid->vk * cid->vk / acos(0.0); double cs0 = 0.25 * cid->vk * cid->vk * cid->vk / acos(0.0);
double csch = 2.0 * cid->vk * cid->sqsfi / cid->c1->gcs; double csch = 2.0 * cid->vk * cid->sqsfi / cid->c1->gcs;
...@@ -1642,6 +1622,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo ...@@ -1642,6 +1622,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
nvtxRangePush("Angle loop"); nvtxRangePush("Angle loop");
#endif #endif
double th = sa->th; double th = sa->th;
int done_dirs = 0;
for (int jth486 = 1; jth486 <= sa->nth; jth486++) { // OpenMP portable? for (int jth486 = 1; jth486 <= sa->nth; jth486++) { // OpenMP portable?
double ph = sa->ph; double ph = sa->ph;
double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0; double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0;
...@@ -1767,8 +1748,6 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo ...@@ -1767,8 +1748,6 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
vtppoanp->append_line(VirtualBinaryLine(value)); vtppoanp->append_line(VirtualBinaryLine(value));
} }
} }
sprintf(virtual_line, " ENSEMBLE AVERAGE, MODE%2d\n", iavm);
output->append_line(virtual_line);
int jlr = 2; int jlr = 2;
for (int ilr210 = 1; ilr210 <= 2; ilr210++) { for (int ilr210 = 1; ilr210 <= 2; ilr210++) {
int ipol = (ilr210 % 2 == 0) ? 1 : -1; int ipol = (ilr210 % 2 == 0) ? 1 : -1;
...@@ -1784,68 +1763,72 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo ...@@ -1784,68 +1763,72 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
double qschum = imag(s0m) * csch; double qschum = imag(s0m) * csch;
double pschum = real(s0m) * csch; double pschum = real(s0m) * csch;
double s0magm = cabs(s0m) * cs0; double s0magm = cabs(s0m) * cs0;
if (inpol == 0) { // if (inpol == 0) {
sprintf(virtual_line, " LIN %2d\n", ipol); // sprintf(virtual_line, " LIN %2d\n", ipol);
output->append_line(virtual_line); // output->append_line(virtual_line);
} else { // label 206 // } else { // label 206
sprintf(virtual_line, " CIRC %2d\n", ipol); // sprintf(virtual_line, " CIRC %2d\n", ipol);
output->append_line(virtual_line); // output->append_line(virtual_line);
} // }
// label 208 // label 208
sprintf(virtual_line, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
output->append_line(virtual_line);
sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n", scasm, abssm, extsm, albdm);
output->append_line(virtual_line);
double alamb = 2.0 * 3.141592653589793238 / cid->vk;
sprintf(virtual_line, "INSERTION: SCASECM %5d%15.7E%15.7E%15.7E%15.7E\n", ipol, alamb, scasm, abssm, extsm);
output->append_line(virtual_line);
sprintf(virtual_line, " ---- SCS/GS -- ABC/GS -- EXS/GS ---\n");
output->append_line(virtual_line);
sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n", qscam, qabsm, qextm);
output->append_line(virtual_line);
sprintf(
virtual_line, " FSAS(%1d,%1d)=%15.7lE%15.7lE FSAS(%1d,%1d)=%15.7lE%15.7lE\n",
ilr210, ilr210, real(cid->c1->fsacm[ilr210 - 1][ilr210 - 1]),
imag(cid->c1->fsacm[ilr210 - 1][ilr210 - 1]), jlr, ilr210,
real(cid->c1->fsacm[jlr - 1][ilr210 - 1]), imag(cid->c1->fsacm[jlr - 1][ilr210 - 1])
);
output->append_line(virtual_line);
sprintf(virtual_line, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschum, pschum, s0magm);
output->append_line(virtual_line);
double rapr = cid->c1->ecscm[ilr210 - 1] - cid->gapm[2][ilr210 - 1]; double rapr = cid->c1->ecscm[ilr210 - 1] - cid->gapm[2][ilr210 - 1];
double cosav = cid->gapm[2][ilr210 - 1] / cid->c1->scscm[ilr210 - 1]; double cosav = cid->gapm[2][ilr210 - 1] / cid->c1->scscm[ilr210 - 1];
double fz = rapr; double fz = rapr;
sprintf(virtual_line, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr); if (ipol == -1) {
output->append_line(virtual_line); output->vec_scs1[jindex - 1] = scasm;
sprintf(virtual_line, " Fk=%15.7lE\n", fz); output->vec_abs1[jindex - 1] = abssm;
output->append_line(virtual_line); output->vec_exs1[jindex - 1] = extsm;
output->vec_albeds1[jindex - 1] = albdm;
output->vec_scsrt1[jindex - 1] = qscam;
output->vec_absrt1[jindex - 1] = qabsm;
output->vec_exsrt1[jindex - 1] = qextm;
output->vec_fsas11[jindex - 1] = cid->c1->fsacm[0][0];
output->vec_fsas21[jindex - 1] = cid->c1->fsacm[1][0];
output->vec_qschu1[jindex -1] = qschum;
output->vec_pschu1[jindex -1] = pschum;
output->vec_s0mag1[jindex -1] = s0magm;
output->vec_cosav1[jindex -1] = cosav;
output->vec_raprs1[jindex -1] = rapr;
output->vec_fk1[jindex -1] = fz;
} else if (ipol == 1) {
output->vec_scs2[jindex - 1] = scasm;
output->vec_abs2[jindex - 1] = abssm;
output->vec_exs2[jindex - 1] = extsm;
output->vec_albeds2[jindex - 1] = albdm;
output->vec_scsrt2[jindex - 1] = qscam;
output->vec_absrt2[jindex - 1] = qabsm;
output->vec_exsrt2[jindex - 1] = qextm;
output->vec_fsas22[jindex - 1] = cid->c1->fsacm[1][1];
output->vec_fsas12[jindex - 1] = cid->c1->fsacm[0][1];
output->vec_qschu2[jindex -1] = qschum;
output->vec_pschu2[jindex -1] = pschum;
output->vec_s0mag2[jindex -1] = s0magm;
output->vec_cosav2[jindex -1] = cosav;
output->vec_raprs2[jindex -1] = rapr;
output->vec_fk2[jindex -1] = fz;
}
} // ilr210 loop } // ilr210 loop
double rmbrif = (real(cid->c1->fsacm[0][0]) - real(cid->c1->fsacm[1][1])) / real(cid->c1->fsacm[0][0]);
double rmdchr = (imag(cid->c1->fsacm[0][0]) - imag(cid->c1->fsacm[1][1])) / imag(cid->c1->fsacm[0][0]);
sprintf(virtual_line, " (RE(FSAS(1,1))-RE(FSAS(2,2)))/RE(FSAS(1,1))=%15.7lE\n", rmbrif);
output->append_line(virtual_line);
sprintf(virtual_line, " (IM(FSAS(1,1))-IM(FSAS(2,2)))/IM(FSAS(1,1))=%15.7lE\n", rmdchr);
output->append_line(virtual_line);
} }
// label 212 // label 212
sprintf(virtual_line, "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n", jth486, jph484, jths, jphs); output->vec_dir_scand[done_dirs] = cid->scan;
output->append_line(virtual_line); output->vec_dir_cfmp[done_dirs] = cid->cfmp;
sprintf(virtual_line, " TIDG=%10.3lE, PIDG=%10.3lE, TSDG=%10.3lE, PSDG=%10.3lE\n", th, ph, ths, phs); output->vec_dir_sfmp[done_dirs] = cid->sfmp;
output->append_line(virtual_line); output->vec_dir_cfsp[done_dirs] = cid->cfsp;
sprintf(virtual_line, " SCAND=%10.3lE\n", cid->scan); output->vec_dir_sfsp[done_dirs] = cid->sfsp;
output->append_line(virtual_line);
sprintf(virtual_line, " CFMP=%15.7lE, SFMP=%15.7lE\n", cid->cfmp, cid->sfmp);
output->append_line(virtual_line);
sprintf(virtual_line, " CFSP=%15.7lE, SFSP=%15.7lE\n", cid->cfsp, cid->sfsp);
output->append_line(virtual_line);
if (isam >= 0) { if (isam >= 0) {
sprintf(virtual_line, " UNI=(%12.5lE,%12.5lE,%12.5lE)\n", cid->un[0], cid->un[1], cid->un[2]); output->vec_dir_un[3 * done_dirs] = cid->un[0];
output->append_line(virtual_line); output->vec_dir_un[3 * done_dirs + 1] = cid->un[1];
sprintf(virtual_line, " UNS=(%12.5lE,%12.5lE,%12.5lE)\n", cid->uns[0], cid->uns[1], cid->uns[2]); output->vec_dir_un[3 * done_dirs + 2] = cid->un[2];
output->append_line(virtual_line); output->vec_dir_uns[3 * done_dirs] = cid->uns[0];
output->vec_dir_uns[3 * done_dirs + 1] = cid->uns[1];
output->vec_dir_uns[3 * done_dirs + 2] = cid->uns[2];
} else { // label 214 } else { // label 214
sprintf(virtual_line, " UN=(%12.5lE,%12.5lE,%12.5lE)\n\n", cid->un[0], cid->un[1], cid->un[2]); output->vec_dir_un[3 * done_dirs] = cid->un[0];
output->append_line(virtual_line); output->vec_dir_un[3 * done_dirs + 1] = cid->un[1];
output->vec_dir_un[3 * done_dirs + 2] = cid->un[2];
output->vec_dir_uns[3 * done_dirs] = 0.0;
output->vec_dir_uns[3 * done_dirs + 1] = 0.0;
output->vec_dir_uns[3 * done_dirs + 2] = 0.0;
} }
// label 220 // label 220
pcros(cid->vk, exri, cid->c1); pcros(cid->vk, exri, cid->c1);
...@@ -1926,8 +1909,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo ...@@ -1926,8 +1909,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
vtppoanp->append_line(VirtualBinaryLine(value)); vtppoanp->append_line(VirtualBinaryLine(value));
} }
} }
sprintf(virtual_line, " SINGLE SCATTERER\n"); oindex = (jindex - 1) * ndirs + done_dirs;
output->append_line(virtual_line);
int jlr = 2; int jlr = 2;
for (int ilr290 = 1; ilr290 <= 2; ilr290++) { for (int ilr290 = 1; ilr290 <= 2; ilr290++) {
int ipol = (ilr290 % 2 == 0) ? 1 : -1; int ipol = (ilr290 % 2 == 0) ? 1 : -1;
...@@ -1943,53 +1925,38 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo ...@@ -1943,53 +1925,38 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
double qschu = imag(s0) * csch; double qschu = imag(s0) * csch;
double pschu = real(s0) * csch; double pschu = real(s0) * csch;
double s0mag = cabs(s0) * cs0; double s0mag = cabs(s0) * cs0;
if (inpol == 0) {
sprintf(virtual_line, " LIN %2d\n", ipol);
output->append_line(virtual_line);
} else { // label 273
sprintf(virtual_line, " CIRC %2d\n", ipol);
output->append_line(virtual_line);
}
// label 275 // label 275
sprintf(virtual_line, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n"); if (ipol == -1) {
output->append_line(virtual_line); output->vec_dir_scs1[oindex] = scasec;
sprintf( output->vec_dir_abs1[oindex] = abssec;
virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n", output->vec_dir_exs1[oindex] = extsec;
scasec, abssec, extsec, albedc output->vec_dir_albeds1[oindex] = albedc;
); output->vec_dir_scsrt1[oindex] = qsca;
output->append_line(virtual_line); output->vec_dir_absrt1[oindex] = qabs;
double alam = 2.0 * 3.141592653589793238 / cid->vk; output->vec_dir_exsrt1[oindex] = qext;
sprintf(virtual_line, "INSERTION: SCASEC %5d%14.7lE%14.7lE%14.7lE%14.7lE\n", ipol, alam, scasec, abssec, extsec); output->vec_dir_fsas11[oindex] = cid->c1->fsac[0][0];
sprintf(virtual_line, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n"); output->vec_dir_fsas21[oindex] = cid->c1->fsac[1][0];
output->append_line(virtual_line); output->vec_dir_sas11[oindex] = cid->c1->sac[0][0];
sprintf( output->vec_dir_sas21[oindex] = cid->c1->sac[1][0];
virtual_line, " %14.7lE%15.7lE%15.7lE\n", output->vec_dir_qschu1[oindex] = qschu;
qsca, qabs, qext output->vec_dir_pschu1[oindex] = pschu;
); output->vec_dir_s0mag1[oindex] = s0mag;
output->append_line(virtual_line); } else if (ipol == 1) {
sprintf( output->vec_dir_scs2[oindex] = scasec;
virtual_line, output->vec_dir_abs2[oindex] = abssec;
" FSAS(%1d,%1d)=%15.7lE%15.7lE FSAS(%1d,%1d)=%15.7lE%15.7lE\n", output->vec_dir_exs2[oindex] = extsec;
ilr290, ilr290, real(cid->c1->fsac[ilr290 - 1][ilr290 - 1]), output->vec_dir_albeds2[oindex] = albedc;
imag(cid->c1->fsac[ilr290 - 1][ilr290 - 1]), jlr, ilr290, output->vec_dir_scsrt2[oindex] = qsca;
real(cid->c1->fsac[jlr - 1][ilr290 - 1]), output->vec_dir_absrt2[oindex] = qabs;
imag(cid->c1->fsac[jlr - 1][ilr290 - 1]) output->vec_dir_exsrt2[oindex] = qext;
); output->vec_dir_fsas22[oindex] = cid->c1->fsac[1][1];
output->append_line(virtual_line); output->vec_dir_fsas12[oindex] = cid->c1->fsac[0][1];
sprintf( output->vec_dir_sas22[oindex] = cid->c1->sac[1][1];
virtual_line, output->vec_dir_sas12[oindex] = cid->c1->sac[0][1];
" SAS(%1d,%1d)=%15.7lE%15.7lE SAS(%1d,%1d)=%15.7lE%15.7lE\n", output->vec_dir_qschu2[oindex] = qschu;
ilr290, ilr290, real(cid->c1->sac[ilr290 - 1][ilr290 - 1]), output->vec_dir_pschu2[oindex] = pschu;
imag(cid->c1->sac[ilr290 - 1][ilr290 - 1]), jlr, ilr290, output->vec_dir_s0mag2[oindex] = s0mag;
real(cid->c1->sac[jlr - 1][ilr290 - 1]), }
imag(cid->c1->sac[jlr - 1][ilr290 - 1])
);
output->append_line(virtual_line);
sprintf(
virtual_line, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
qschu, pschu, s0mag
);
output->append_line(virtual_line);
bool goto290 = isam >= 0 && (jths > 1 || jphs > 1); bool goto290 = isam >= 0 && (jths > 1 || jphs > 1);
if (!goto290) { if (!goto290) {
cid->gapv[0] = cid->gap[0][ilr290 - 1]; cid->gapv[0] = cid->gap[0][ilr290 - 1];
...@@ -1999,12 +1966,25 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo ...@@ -1999,12 +1966,25 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
double scatts = cid->c1->scsc[ilr290 - 1]; double scatts = cid->c1->scsc[ilr290 - 1];
double rapr, cosav, fp, fn, fk, fx, fy, fz; double rapr, cosav, fp, fn, fk, fx, fy, fz;
rftr(cid->u, cid->up, cid->un, cid->gapv, extins, scatts, rapr, cosav, fp, fn, fk, fx, fy, fz); rftr(cid->u, cid->up, cid->un, cid->gapv, extins, scatts, rapr, cosav, fp, fn, fk, fx, fy, fz);
sprintf(virtual_line, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr); if (ipol == -1) {
output->append_line(virtual_line); output->vec_dir_cosav1[oindex] = cosav;
sprintf(virtual_line, " Fl=%15.7lE, Fr=%15.7lE, Fk=%15.7lE\n", fp, fn, fk); output->vec_dir_rapr1[oindex] = rapr;
output->append_line(virtual_line); output->vec_dir_fl1[oindex] = fp;
sprintf(virtual_line, " Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n", fx, fy, fz); output->vec_dir_fr1[oindex] = fn;
output->append_line(virtual_line); output->vec_dir_fk1[oindex] = fk;
output->vec_dir_fx1[oindex] = fx;
output->vec_dir_fy1[oindex] = fy;
output->vec_dir_fz1[oindex] = fz;
} else if (ipol == 1) {
output->vec_dir_cosav2[oindex] = cosav;
output->vec_dir_rapr2[oindex] = rapr;
output->vec_dir_fl2[oindex] = fp;
output->vec_dir_fr2[oindex] = fn;
output->vec_dir_fk2[oindex] = fk;
output->vec_dir_fx2[oindex] = fx;
output->vec_dir_fy2[oindex] = fy;
output->vec_dir_fz2[oindex] = fz;
}
cid->tqev[0] = cid->tqce[ilr290 - 1][0]; cid->tqev[0] = cid->tqce[ilr290 - 1][0];
cid->tqev[1] = cid->tqce[ilr290 - 1][1]; cid->tqev[1] = cid->tqce[ilr290 - 1][1];
cid->tqev[2] = cid->tqce[ilr290 - 1][2]; cid->tqev[2] = cid->tqce[ilr290 - 1][2];
...@@ -2013,45 +1993,48 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo ...@@ -2013,45 +1993,48 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
cid->tqsv[2] = cid->tqcs[ilr290 - 1][2]; cid->tqsv[2] = cid->tqcs[ilr290 - 1][2];
double tep, ten, tek, tsp, tsn, tsk; double tep, ten, tek, tsp, tsn, tsk;
tqr(cid->u, cid->up, cid->un, cid->tqev, cid->tqsv, tep, ten, tek, tsp, tsn, tsk); tqr(cid->u, cid->up, cid->un, cid->tqev, cid->tqsv, tep, ten, tek, tsp, tsn, tsk);
sprintf(virtual_line, " TQEl=%15.7lE, TQEr=%15.7lE, TQEk=%15.7lE\n", tep, ten, tek); if (ipol == -1) {
output->append_line(virtual_line); output->vec_dir_tqel1[oindex] = tep;
sprintf(virtual_line, " TQSl=%15.7lE, TQSr=%15.7lE, TQSk=%15.7lE\n", tsp, tsn, tsk); output->vec_dir_tqer1[oindex] = ten;
output->append_line(virtual_line); output->vec_dir_tqek1[oindex] = tek;
sprintf( output->vec_dir_tqsl1[oindex] = tsp;
virtual_line, " TQEx=%15.7lE, TQEy=%15.7lE, TQEz=%15.7lE\n", output->vec_dir_tqsr1[oindex] = tsn;
cid->tqce[ilr290 - 1][0], cid->tqce[ilr290 - 1][1], cid->tqce[ilr290 - 1][2] output->vec_dir_tqsk1[oindex] = tsk;
); output->vec_dir_tqex1[oindex] = cid->tqce[0][0];
output->append_line(virtual_line); output->vec_dir_tqey1[oindex] = cid->tqce[0][1];
sprintf( output->vec_dir_tqez1[oindex] = cid->tqce[0][2];
virtual_line, " TQSx=%15.7lE, TQSy=%15.7lE, TQSz=%15.7lE\n", output->vec_dir_tqsx1[oindex] = cid->tqcs[0][0];
cid->tqcs[ilr290 - 1][0], cid->tqcs[ilr290 - 1][1], cid->tqcs[ilr290 - 1][2] output->vec_dir_tqsy1[oindex] = cid->tqcs[0][1];
); output->vec_dir_tqsz1[oindex] = cid->tqcs[0][2];
output->append_line(virtual_line); } else if (ipol == 1) {
output->vec_dir_tqel2[oindex] = tep;
output->vec_dir_tqer2[oindex] = ten;
output->vec_dir_tqek2[oindex] = tek;
output->vec_dir_tqsl2[oindex] = tsp;
output->vec_dir_tqsr2[oindex] = tsn;
output->vec_dir_tqsk2[oindex] = tsk;
output->vec_dir_tqex2[oindex] = cid->tqce[1][0];
output->vec_dir_tqey2[oindex] = cid->tqce[1][1];
output->vec_dir_tqez2[oindex] = cid->tqce[1][2];
output->vec_dir_tqsx2[oindex] = cid->tqcs[1][0];
output->vec_dir_tqsy2[oindex] = cid->tqcs[1][1];
output->vec_dir_tqsz2[oindex] = cid->tqcs[1][2];
}
} }
} //ilr290 loop } //ilr290 loop
double rbirif = (real(cid->c1->fsac[0][0]) - real(cid->c1->fsac[1][1])) / real(cid->c1->fsac[0][0]);
double rdichr = (imag(cid->c1->fsac[0][0]) - imag(cid->c1->fsac[1][1])) / imag(cid->c1->fsac[0][0]);
sprintf(virtual_line, " (RE(FSAS(1,1))-RE(FSAS(2,2)))/RE(FSAS(1,1))=%15.7lE\n", rbirif);
output->append_line(virtual_line);
sprintf(virtual_line, " (IM(FSAS(1,1))-IM(FSAS(2,2)))/IM(FSAS(1,1))=%15.7lE\n", rdichr);
output->append_line(virtual_line);
sprintf(virtual_line, " MULL\n");
output->append_line(virtual_line);
for (int i = 0; i < 4; i++) { for (int i = 0; i < 4; i++) {
sprintf( oindex = 16 * (jindex - 1) * ndirs + 16 * done_dirs + 4 * i;
virtual_line, " %15.7lE%15.7lE%15.7lE%15.7lE\n", output->vec_dir_mull[oindex] = cid->cmul[i][0];
cid->cmul[i][0], cid->cmul[i][1], cid->cmul[i][2], cid->cmul[i][3] output->vec_dir_mull[oindex + 1] = cid->cmul[i][1];
); output->vec_dir_mull[oindex + 2] = cid->cmul[i][2];
output->append_line(virtual_line); output->vec_dir_mull[oindex + 3] = cid->cmul[i][3];
} }
sprintf(virtual_line, " MULLLR\n");
output->append_line(virtual_line);
for (int i = 0; i < 4; i++) { for (int i = 0; i < 4; i++) {
sprintf( oindex = 16 * (jindex - 1) * ndirs + 16 * done_dirs + 4 * i;
virtual_line, " %15.7lE%15.7lE%15.7lE%15.7lE\n", output->vec_dir_mulllr[oindex] = cid->cmullr[i][0];
cid->cmullr[i][0], cid->cmullr[i][1], cid->cmullr[i][2], cid->cmullr[i][3] output->vec_dir_mulllr[oindex + 1] = cid->cmullr[i][1];
); output->vec_dir_mulllr[oindex + 2] = cid->cmullr[i][2];
output->append_line(virtual_line); output->vec_dir_mulllr[oindex + 3] = cid->cmullr[i][3];
} }
if (iavm != 0) { if (iavm != 0) {
mmulc(cid->c1->vintm, cid->cmullr, cid->cmul); mmulc(cid->c1->vintm, cid->cmullr, cid->cmul);
...@@ -2069,37 +2052,25 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo ...@@ -2069,37 +2052,25 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
vtppoanp->append_line(VirtualBinaryLine(value)); vtppoanp->append_line(VirtualBinaryLine(value));
} }
} }
sprintf(virtual_line, " ENSEMBLE AVERAGE, MODE%2d\n", iavm);
output->append_line(virtual_line);
if (inpol == 0) {
sprintf(virtual_line, " LIN\n");
output->append_line(virtual_line);
} else { // label 316
sprintf(virtual_line, " CIRC\n");
output->append_line(virtual_line);
}
// label 318 // label 318
sprintf(virtual_line, " MULL\n");
output->append_line(virtual_line);
for (int i = 0; i < 4; i++) { for (int i = 0; i < 4; i++) {
sprintf( oindex = 16 * (jindex - 1) + 4 * i; // if IAVM fails, try adding directions
virtual_line, " %15.7lE%15.7lE%15.7lE%15.7lE\n", output->vec_dir_mull[oindex] = cid->cmul[i][0];
cid->cmul[i][0], cid->cmul[i][1], cid->cmul[i][2], cid->cmul[i][3] output->vec_dir_mull[oindex + 1] = cid->cmul[i][1];
); output->vec_dir_mull[oindex + 2] = cid->cmul[i][2];
output->append_line(virtual_line); output->vec_dir_mull[oindex + 3] = cid->cmul[i][3];
} }
sprintf(virtual_line, " MULLLR\n");
output->append_line(virtual_line);
for (int i = 0; i < 4; i++) { for (int i = 0; i < 4; i++) {
sprintf( oindex = 16 * (jindex - 1) + 4 * i; // if IAVM fails, try adding directions
virtual_line, " %15.7lE%15.7lE%15.7lE%15.7lE\n", output->vec_dir_mulllr[oindex] = cid->cmul[i][0];
cid->cmullr[i][0], cid->cmullr[i][1], cid->cmullr[i][2], cid->cmullr[i][3] output->vec_dir_mulllr[oindex + 1] = cid->cmullr[i][1];
); output->vec_dir_mulllr[oindex + 2] = cid->cmullr[i][2];
output->append_line(virtual_line); output->vec_dir_mulllr[oindex + 3] = cid->cmullr[i][3];
} }
} }
// label 420, continues jphs loop // label 420, continues jphs loop
if (isam < 1) phs += sa->phsstp; if (isam < 1) phs += sa->phsstp;
done_dirs++;
} // jphs loop, labeled 480 } // jphs loop, labeled 480
if (isam <= 1) thsl += sa->thsstp; if (isam <= 1) thsl += sa->thsstp;
} // jths loop, labeled 482 } // jths loop, labeled 482
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment