From 1c5cd21511ca988def355105c0e4ba01cc6a37a0 Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 24 Jul 2025 16:31:34 +0200 Subject: [PATCH 1/7] Update Phonon_simple fixing and modernising Make sure that vf_list has length of 100 and allow the user to choose the amount of scan steps actually taken, if they need to optimise the instrument. This counteracts overflow of vf_list. Also change from parms array to parms struct in order to improve readability. --- mcstas-comps/samples/Phonon_simple.comp | 165 +++++++++++++++--------- 1 file changed, 103 insertions(+), 62 deletions(-) diff --git a/mcstas-comps/samples/Phonon_simple.comp b/mcstas-comps/samples/Phonon_simple.comp index 7e5ea0fec8..ad47cc1a25 100644 --- a/mcstas-comps/samples/Phonon_simple.comp +++ b/mcstas-comps/samples/Phonon_simple.comp @@ -66,7 +66,9 @@ DEFINE COMPONENT Phonon_simple SETTING PARAMETERS (radius,yheight,sigma_abs,sigma_inc,a,b,M,c,DW,T, -target_x=0, target_y=0, target_z=0, int target_index=0,focus_r=0,focus_xw=0,focus_yh=0,focus_aw=0,focus_ah=0, gap=0) + target_x=0, target_y=0, target_z=0, int target_index=0, + focus_r=0,focus_xw=0,focus_yh=0,focus_aw=0,focus_ah=0, + gap=0, int e_steps_low=50, int e_steps_high=50) /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ SHARE @@ -75,6 +77,24 @@ SHARE #define PHONON_SIMPLE $Revision$ #define T2E (1/11.605) /* Kelvin to meV */ + struct neutron_and_phonon_params { + // Statically allocate vectors that are always 3 + double vf; // Final velocity size + double vi; // Initial velocity size + double vv_x; // vv is the unit vector of the final velocity vector + double vv_y; + double vv_z; + double vi_x; // vi is the initial velocity vector + double vi_y; + double vi_z; + double a_; // d spacing of the cubic lattice + double c_; // Speed of sound in the material + double gap_; // optional spin gap + double ah; // Half of a + int e_steps_high_; + int e_steps_low_; + }; + #pragma acc routine double nbose(double omega, double T) /* Other name ?? */ { @@ -103,7 +123,7 @@ void fatalerror(char *s) } #pragma acc routine - double omega_q(double* parms) + double omega_q(struct neutron_and_phonon_params *parms) { /* dispersion in units of meV */ double vi, vf, vv_x, vv_y, vv_z, vi_x, vi_y, vi_z; @@ -111,18 +131,18 @@ void fatalerror(char *s) double ah, a, c; double gap; - vf=parms[0]; - vi=parms[1]; - vv_x=parms[2]; - vv_y=parms[3]; - vv_z=parms[4]; - vi_x=parms[5]; - vi_y=parms[6]; - vi_z=parms[7]; - a =parms[8]; - c =parms[9]; - gap =parms[10]; - ah=a/2.0; + vf=parms->vf; + vi=parms->vi; + vv_x=parms->vv_x; + vv_y=parms->vv_y; + vv_z=parms->vv_z; + vi_x=parms->vi_x; + vi_y=parms->vi_y; + vi_z=parms->vi_z; + a =parms->a_; + c =parms->c_; + gap =parms->gap_; + ah =parms->ah; qx=V2K*(vi_x-vf*vv_x); qy=V2K*(vi_y-vf*vv_y); @@ -141,15 +161,15 @@ void fatalerror(char *s) } -double zridd(double (*func)(double*), double x1, double x2, double *parms, double xacc) +double zridd(double (*func)(struct neutron_and_phonon_params*), double x1, double x2, struct neutron_and_phonon_params *parms, double xacc) { int j; double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; - parms[0]=x1; - fl=(*func)(parms); - parms[0]=x2; - fh=(*func)(parms); + parms->vf=x1; + fl=func(parms); + parms->vf=x2; + fh=func(parms); if (fl*fh >= 0) { if (fl==0) return x1; @@ -164,8 +184,8 @@ double zridd(double (*func)(double*), double x1, double x2, double *parms, doubl for (j=1; jvf=xm; + fm=func(parms); s=sqrt(fm*fm-fl*fh); if (s == 0.0) return ans; @@ -173,8 +193,8 @@ double zridd(double (*func)(double*), double x1, double x2, double *parms, doubl if (fabs(xnew-ans) <= xacc) return ans; ans=xnew; - parms[0]=ans; - fnew=(*func)(parms); + parms->vf=ans; + fnew=func(parms); if (fnew == 0.0) return ans; if (fabs(fm)*SIGN(fnew) != fm) { @@ -206,14 +226,14 @@ double zridd(double (*func)(double*), double x1, double x2, double *parms, doubl } #pragma acc routine -double zridd_gpu(double x1, double x2, double *parms, double xacc) +double zridd_gpu(double x1, double x2, struct neutron_and_phonon_params *parms, double xacc) { int j; double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; - parms[0]=x1; + parms->vf=x1; fl=omega_q(parms); - parms[0]=x2; + parms->vf=x2; fh=omega_q(parms); if (fl*fh >= 0) { @@ -229,7 +249,7 @@ double zridd_gpu(double x1, double x2, double *parms, double xacc) for (j=1; jvf=xm; fm=omega_q(parms); s=sqrt(fm*fm-fl*fh); if (s == 0.0) @@ -238,7 +258,7 @@ double zridd_gpu(double x1, double x2, double *parms, double xacc) if (fabs(xnew-ans) <= xacc) return ans; ans=xnew; - parms[0]=ans; + parms->vf=ans; fnew=omega_q(parms); if (fnew == 0.0) return ans; if (fabs(fm)*SIGN(fnew) != fm) @@ -273,30 +293,40 @@ double zridd_gpu(double x1, double x2, double *parms, double xacc) #define ROOTACC 1e-8 - int findroots(double brack_low, double brack_mid, double brack_high, double *list, int* index, double (*f)(double*), double *parms) + void findroots(double brack_low, double brack_mid, double brack_high, + double *list, int* index, double (*f)(struct neutron_and_phonon_params*), + struct neutron_and_phonon_params *parms) { - double root,range=brack_mid-brack_low; - int i, steps=100; - - for (i=0; ie_steps_low_; i++){ + root = zridd(f, brack_low+range_low*i/ parms->e_steps_low_, + brack_low+range_low*(i+1)/ parms->e_steps_low_, + parms, ROOTACC); if (root != UNUSED) { list[(*index)++]=root; } } - root = zridd(f, brack_mid, brack_high, (double *)parms, ROOTACC); - if (root != UNUSED) - { + + // Then in energy gain for the neutron + for (int i=0; ie_steps_high_; i++){ + root = zridd(f, brack_low+range_high*i/ parms->e_steps_high_, + brack_low+range_high*(i+1)/ parms->e_steps_high_, + parms, ROOTACC); + if (root != UNUSED){ list[(*index)++]=root; } + } } #pragma acc routine - int findroots_gpu(double brack_low, double brack_mid, double brack_high, double *list, int* index, double *parms) + void findroots_gpu(double brack_low, double brack_mid, double brack_high, + double *list, int* index, struct neutron_and_phonon_params *parms) { double root,range=brack_mid-brack_low; int i, steps=100; @@ -305,13 +335,13 @@ double zridd_gpu(double x1, double x2, double *parms, double xacc) { root = zridd_gpu(brack_low+range*i/(int)steps, brack_low+range*(i+1)/(int)steps, - (double *)parms, ROOTACC); + parms, ROOTACC); if (root != UNUSED) { list[(*index)++]=root; } } - root = zridd_gpu(brack_mid, brack_high, (double *)parms, ROOTACC); + root = zridd_gpu(brack_mid, brack_high, parms, ROOTACC); if (root != UNUSED) { list[(*index)++]=root; @@ -329,13 +359,25 @@ DECLARE double V_my_s; double V_my_a_v; double DV; + double vf_list[100];// List of allowed final velocities. Has length of scan_steps + struct neutron_and_phonon_params parms; + %} INITIALIZE %{ + memset(vf_list, 0, sizeof(vf_list)); V_rho = 4/(a*a*a); V_my_s = (V_rho * 100 * sigma_inc); V_my_a_v = (V_rho * 100 * sigma_abs * 2200); DV = 0.001; /* Velocity change used for numerical derivative */ + + // Set constant parameters for parms object + parms.a_ = a; + parms.c_ = c; + parms.gap_ = gap; + parms.ah = a/2.0; + parms.e_steps_high_ = e_steps_high; + parms.e_steps_low_ = e_steps_low; /* now compute target coords if a component index is supplied */ if (!target_index && !target_x && !target_y && !target_z) target_index=1; @@ -368,11 +410,9 @@ TRACE double bose_factor; /* Calculated value of the Bose factor */ double omega; /* energy transfer */ int nf, index; /* Number of allowed final velocities */ - double vf_list[2]; /* List of allowed final velocities */ double J_factor; /* Jacobian from delta fnc.s in cross section */ double f1, f2; /* probed values of omega_q minus omega */ double p1,p2,p3,p4,p5; /* temporary multipliers */ - double parms[11]; if(cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius, yheight)) { @@ -405,29 +445,27 @@ TRACE } NORM(vx, vy, vz); nf=0; - parms[0]=-1; - parms[1]=v_i; - parms[2]=vx; - parms[3]=vy; - parms[4]=vz; - parms[5]=vx_i; - parms[6]=vy_i; - parms[7]=vz_i; - parms[8]=a; - parms[9]=c; - parms[10]=gap; + parms.vf = -1; + parms.vi = v_i; + parms.vv_x = vx; + parms.vv_y = vy; + parms.vv_z = vz; + parms.vi_x = vx_i; + parms.vi_y = vy_i; + parms.vi_z = vz_i; + #ifndef OPENACC - findroots(0, v_i, v_i+2*c*V2K/VS2E, vf_list, &nf, omega_q, parms); + findroots(0, v_i, v_i+2*c*V2K/VS2E, vf_list, &nf, omega_q, &parms); #else - findroots_gpu(0, v_i, v_i+2*c*V2K/VS2E, vf_list, &nf, parms); + findroots_gpu(0, v_i, v_i+2*c*V2K/VS2E, vf_list, &nf, &parms); #endif index=(int)floor(rand01()*nf); if (index<2) { v_f=vf_list[index]; - parms[0]=v_f-DV; - f1=omega_q(parms); - parms[0]=v_f+DV; - f2=omega_q(parms); + parms.vf=v_f-DV; + f1=omega_q(&parms); + parms.vf=v_f+DV; + f2=omega_q(&parms); J_factor = fabs(f2-f1)/(2*DV); omega=VS2E*(v_i*v_i-v_f*v_f); vx *= v_f; @@ -460,6 +498,9 @@ TRACE ABSORB; // findroots returned junk } } /* else transmit: Neutron did not hit the sample */ + + // Finally, set vf_list to zero again + memset(vf_list, 0, nf * sizeof(double)); %} MCDISPLAY From 4f1f038a79a431762b3d68c5e33c1717cbdb9bfe Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 14 Aug 2025 10:11:42 +0200 Subject: [PATCH 2/7] Change vf_list to be as long as the energy scans --- mcstas-comps/samples/Phonon_simple.comp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/mcstas-comps/samples/Phonon_simple.comp b/mcstas-comps/samples/Phonon_simple.comp index ad47cc1a25..c24ad34abd 100644 --- a/mcstas-comps/samples/Phonon_simple.comp +++ b/mcstas-comps/samples/Phonon_simple.comp @@ -51,6 +51,9 @@ * target_z: [m] position of target to focus at. Straight ahead. * target_index: [1] relative index of component to focus at, e.g. next is +1 * gap: [meV] Bandgap energy (unphysical) +* e_steps_low [50] Amount of possible intersections beneath the elastic line +* e_steps_high [50] Amount of possible intersections above the elastic line +* * * CALCULATED PARAMETERS: * V_rho: [AA^-3] Atomic density @@ -359,13 +362,13 @@ DECLARE double V_my_s; double V_my_a_v; double DV; - double vf_list[100];// List of allowed final velocities. Has length of scan_steps + double* vf_list;// List of allowed final velocities. Has length of scan_steps struct neutron_and_phonon_params parms; %} INITIALIZE %{ - memset(vf_list, 0, sizeof(vf_list)); + vf_list = (double *)calloc( e_steps_low*e_steps_high, sizeof(double)); V_rho = 4/(a*a*a); V_my_s = (V_rho * 100 * sigma_inc); V_my_a_v = (V_rho * 100 * sigma_abs * 2200); From 27901e2c6bbf90b7c9c1420c7cf4d125038524c8 Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 14 Aug 2025 11:03:30 +0200 Subject: [PATCH 3/7] Remove limitation that Phonon_simple only chooses first solution --- mcstas-comps/samples/Phonon_simple.comp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/mcstas-comps/samples/Phonon_simple.comp b/mcstas-comps/samples/Phonon_simple.comp index c24ad34abd..8e6eeda200 100644 --- a/mcstas-comps/samples/Phonon_simple.comp +++ b/mcstas-comps/samples/Phonon_simple.comp @@ -463,7 +463,7 @@ TRACE findroots_gpu(0, v_i, v_i+2*c*V2K/VS2E, vf_list, &nf, &parms); #endif index=(int)floor(rand01()*nf); - if (index<2) { + v_f=vf_list[index]; parms.vf=v_f-DV; f1=omega_q(&parms); @@ -497,9 +497,7 @@ TRACE p4 = 2*VS2E*v_f/J_factor; /* Jacobian of delta functions in cross section */ p5 = b*b/M; /* Cross section factor 2 */ p *= p1*p2*p3*p4*p5; - } else { - ABSORB; // findroots returned junk - } + } /* else transmit: Neutron did not hit the sample */ // Finally, set vf_list to zero again From 956c950253c73d3e8f67009c07a52fd807750d75 Mon Sep 17 00:00:00 2001 From: Diablo Date: Mon, 18 Aug 2025 11:45:41 +0200 Subject: [PATCH 4/7] Change from one struct, to two --- mcstas-comps/samples/Phonon_simple.comp | 161 +++++++++++++----------- 1 file changed, 86 insertions(+), 75 deletions(-) diff --git a/mcstas-comps/samples/Phonon_simple.comp b/mcstas-comps/samples/Phonon_simple.comp index 8e6eeda200..213fd10ddc 100644 --- a/mcstas-comps/samples/Phonon_simple.comp +++ b/mcstas-comps/samples/Phonon_simple.comp @@ -80,7 +80,15 @@ SHARE #define PHONON_SIMPLE $Revision$ #define T2E (1/11.605) /* Kelvin to meV */ - struct neutron_and_phonon_params { + struct phonon_params { + double a_; // d spacing of the cubic lattice + double c_; // Speed of sound in the material + double gap_; // optional spin gap + double ah; // Half of a + int e_steps_high_; + int e_steps_low_; + }; + struct neutron_params { // Statically allocate vectors that are always 3 double vf; // Final velocity size double vi; // Initial velocity size @@ -90,12 +98,6 @@ SHARE double vi_x; // vi is the initial velocity vector double vi_y; double vi_z; - double a_; // d spacing of the cubic lattice - double c_; // Speed of sound in the material - double gap_; // optional spin gap - double ah; // Half of a - int e_steps_high_; - int e_steps_low_; }; #pragma acc routine @@ -126,7 +128,7 @@ void fatalerror(char *s) } #pragma acc routine - double omega_q(struct neutron_and_phonon_params *parms) + double omega_q(struct neutron_params *neutron, struct phonon_params *phonon) { /* dispersion in units of meV */ double vi, vf, vv_x, vv_y, vv_z, vi_x, vi_y, vi_z; @@ -134,18 +136,18 @@ void fatalerror(char *s) double ah, a, c; double gap; - vf=parms->vf; - vi=parms->vi; - vv_x=parms->vv_x; - vv_y=parms->vv_y; - vv_z=parms->vv_z; - vi_x=parms->vi_x; - vi_y=parms->vi_y; - vi_z=parms->vi_z; - a =parms->a_; - c =parms->c_; - gap =parms->gap_; - ah =parms->ah; + vf=neutron->vf; + vi=neutron->vi; + vv_x=neutron->vv_x; + vv_y=neutron->vv_y; + vv_z=neutron->vv_z; + vi_x=neutron->vi_x; + vi_y=neutron->vi_y; + vi_z=neutron->vi_z; + a =phonon->a_; + c =phonon->c_; + gap =phonon->gap_; + ah =phonon->ah; qx=V2K*(vi_x-vf*vv_x); qy=V2K*(vi_y-vf*vv_y); @@ -164,15 +166,17 @@ void fatalerror(char *s) } -double zridd(double (*func)(struct neutron_and_phonon_params*), double x1, double x2, struct neutron_and_phonon_params *parms, double xacc) +double zridd(double (*func)(struct neutron_params*, struct phonon_params*), + double x1, double x2, struct neutron_params *neutron, + struct phonon_params *phonon, double xacc) { int j; double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; - parms->vf=x1; - fl=func(parms); - parms->vf=x2; - fh=func(parms); + neutron->vf=x1; + fl=func(neutron, phonon); + neutron->vf=x2; + fh=func(neutron, phonon); if (fl*fh >= 0) { if (fl==0) return x1; @@ -187,8 +191,8 @@ double zridd(double (*func)(struct neutron_and_phonon_params*), double x1, doubl for (j=1; jvf=xm; - fm=func(parms); + neutron->vf=xm; + fm=func(neutron,phonon); s=sqrt(fm*fm-fl*fh); if (s == 0.0) return ans; @@ -196,8 +200,8 @@ double zridd(double (*func)(struct neutron_and_phonon_params*), double x1, doubl if (fabs(xnew-ans) <= xacc) return ans; ans=xnew; - parms->vf=ans; - fnew=func(parms); + neutron->vf=ans; + fnew=func(neutron,phonon); if (fnew == 0.0) return ans; if (fabs(fm)*SIGN(fnew) != fm) { @@ -229,15 +233,15 @@ double zridd(double (*func)(struct neutron_and_phonon_params*), double x1, doubl } #pragma acc routine -double zridd_gpu(double x1, double x2, struct neutron_and_phonon_params *parms, double xacc) +double zridd_gpu(double x1, double x2, struct neutron_params *neutron, struct phonon_params* phonon, double xacc) { int j; double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; - parms->vf=x1; - fl=omega_q(parms); - parms->vf=x2; - fh=omega_q(parms); + neutron->vf=x1; + fl=omega_q(neutron, phonon); + neutron->vf=x2; + fh=omega_q(neutron, phonon); if (fl*fh >= 0) { if (fl==0) return x1; @@ -252,8 +256,8 @@ double zridd_gpu(double x1, double x2, struct neutron_and_phonon_params *parms, for (j=1; jvf=xm; - fm=omega_q(parms); + neutron->vf=xm; + fm=omega_q(neutron, phonon); s=sqrt(fm*fm-fl*fh); if (s == 0.0) return ans; @@ -261,8 +265,8 @@ double zridd_gpu(double x1, double x2, struct neutron_and_phonon_params *parms, if (fabs(xnew-ans) <= xacc) return ans; ans=xnew; - parms->vf=ans; - fnew=omega_q(parms); + neutron->vf=ans; + fnew=omega_q(neutron, phonon); if (fnew == 0.0) return ans; if (fabs(fm)*SIGN(fnew) != fm) { @@ -297,8 +301,8 @@ double zridd_gpu(double x1, double x2, struct neutron_and_phonon_params *parms, #define ROOTACC 1e-8 void findroots(double brack_low, double brack_mid, double brack_high, - double *list, int* index, double (*f)(struct neutron_and_phonon_params*), - struct neutron_and_phonon_params *parms) + double *list, int* index, double (*f)(struct neutron_params*, struct phonon_params*), + struct neutron_params *neutron, struct phonon_params *phonon) { double root; // Energy gain and energy loss spaces are not equally big. We check uniformly @@ -306,10 +310,10 @@ double zridd_gpu(double x1, double x2, struct neutron_and_phonon_params *parms, double range_low = brack_mid - brack_low; double range_high = brack_high - brack_mid; // First in energy loss for the neutron - for (int i=0; ie_steps_low_; i++){ - root = zridd(f, brack_low+range_low*i/ parms->e_steps_low_, - brack_low+range_low*(i+1)/ parms->e_steps_low_, - parms, ROOTACC); + for (int i=0; ie_steps_low_; i++){ + root = zridd(f, brack_low+range_low*i/ phonon->e_steps_low_, + brack_low+range_low*(i+1)/ phonon->e_steps_low_, + neutron, phonon, ROOTACC); if (root != UNUSED) { list[(*index)++]=root; @@ -317,10 +321,10 @@ double zridd_gpu(double x1, double x2, struct neutron_and_phonon_params *parms, } // Then in energy gain for the neutron - for (int i=0; ie_steps_high_; i++){ - root = zridd(f, brack_low+range_high*i/ parms->e_steps_high_, - brack_low+range_high*(i+1)/ parms->e_steps_high_, - parms, ROOTACC); + for (int i=0; ie_steps_high_; i++){ + root = zridd(f, brack_low+range_high*i/ phonon->e_steps_high_, + brack_low+range_high*(i+1)/ phonon->e_steps_high_, + neutron, phonon, ROOTACC); if (root != UNUSED){ list[(*index)++]=root; } @@ -329,7 +333,7 @@ double zridd_gpu(double x1, double x2, struct neutron_and_phonon_params *parms, #pragma acc routine void findroots_gpu(double brack_low, double brack_mid, double brack_high, - double *list, int* index, struct neutron_and_phonon_params *parms) + double *list, int* index, struct neutron_params *neutron, struct phonon_params *phonon) { double root,range=brack_mid-brack_low; int i, steps=100; @@ -338,13 +342,13 @@ double zridd_gpu(double x1, double x2, struct neutron_and_phonon_params *parms, { root = zridd_gpu(brack_low+range*i/(int)steps, brack_low+range*(i+1)/(int)steps, - parms, ROOTACC); + neutron,phonon, ROOTACC); if (root != UNUSED) { list[(*index)++]=root; } } - root = zridd_gpu(brack_mid, brack_high, parms, ROOTACC); + root = zridd_gpu(brack_mid, brack_high, neutron,phonon, ROOTACC); if (root != UNUSED) { list[(*index)++]=root; @@ -363,7 +367,8 @@ DECLARE double V_my_a_v; double DV; double* vf_list;// List of allowed final velocities. Has length of scan_steps - struct neutron_and_phonon_params parms; + struct neutron_params neutron; + struct phonon_params phonon; %} INITIALIZE @@ -375,12 +380,12 @@ INITIALIZE DV = 0.001; /* Velocity change used for numerical derivative */ // Set constant parameters for parms object - parms.a_ = a; - parms.c_ = c; - parms.gap_ = gap; - parms.ah = a/2.0; - parms.e_steps_high_ = e_steps_high; - parms.e_steps_low_ = e_steps_low; + phonon.a_ = a; + phonon.c_ = c; + phonon.gap_ = gap; + phonon.ah = a/2.0; + phonon.e_steps_high_ = e_steps_high; + phonon.e_steps_low_ = e_steps_low; /* now compute target coords if a component index is supplied */ if (!target_index && !target_x && !target_y && !target_z) target_index=1; @@ -448,27 +453,27 @@ TRACE } NORM(vx, vy, vz); nf=0; - parms.vf = -1; - parms.vi = v_i; - parms.vv_x = vx; - parms.vv_y = vy; - parms.vv_z = vz; - parms.vi_x = vx_i; - parms.vi_y = vy_i; - parms.vi_z = vz_i; + neutron.vf = -1; + neutron.vi = v_i; + neutron.vv_x = vx; + neutron.vv_y = vy; + neutron.vv_z = vz; + neutron.vi_x = vx_i; + neutron.vi_y = vy_i; + neutron.vi_z = vz_i; #ifndef OPENACC - findroots(0, v_i, v_i+2*c*V2K/VS2E, vf_list, &nf, omega_q, &parms); + findroots(0, v_i, v_i+2*c*V2K/VS2E, vf_list, &nf, omega_q, &neutron, &phonon); #else - findroots_gpu(0, v_i, v_i+2*c*V2K/VS2E, vf_list, &nf, &parms); + findroots_gpu(0, v_i, v_i+2*c*V2K/VS2E, vf_list, &nf, &neutron, &phonon); #endif index=(int)floor(rand01()*nf); v_f=vf_list[index]; - parms.vf=v_f-DV; - f1=omega_q(&parms); - parms.vf=v_f+DV; - f2=omega_q(&parms); + neutron.vf=v_f-DV; + f1=omega_q(&neutron, &phonon); + neutron.vf=v_f+DV; + f2=omega_q(&neutron, &phonon); J_factor = fabs(f2-f1)/(2*DV); omega=VS2E*(v_i*v_i-v_f*v_f); vx *= v_f; @@ -497,11 +502,17 @@ TRACE p4 = 2*VS2E*v_f/J_factor; /* Jacobian of delta functions in cross section */ p5 = b*b/M; /* Cross section factor 2 */ p *= p1*p2*p3*p4*p5; - + // Finally, set vf_list to zero again + memset(vf_list, 0, nf * sizeof(double)); } /* else transmit: Neutron did not hit the sample */ - // Finally, set vf_list to zero again - memset(vf_list, 0, nf * sizeof(double)); + +%} + +FINALLY +%{ + free(vf_list); + %} MCDISPLAY From 28d582e7113337868b439d2d5be1d5e898528827 Mon Sep 17 00:00:00 2001 From: Diablo Date: Mon, 18 Aug 2025 13:16:49 +0200 Subject: [PATCH 5/7] Move vf_list and neutron struct to trace in order to comply with GPU and general mcstas logic --- mcstas-comps/samples/Phonon_simple.comp | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/mcstas-comps/samples/Phonon_simple.comp b/mcstas-comps/samples/Phonon_simple.comp index 213fd10ddc..e84fb1f2a7 100644 --- a/mcstas-comps/samples/Phonon_simple.comp +++ b/mcstas-comps/samples/Phonon_simple.comp @@ -366,14 +366,12 @@ DECLARE double V_my_s; double V_my_a_v; double DV; - double* vf_list;// List of allowed final velocities. Has length of scan_steps - struct neutron_params neutron; struct phonon_params phonon; %} INITIALIZE %{ - vf_list = (double *)calloc( e_steps_low*e_steps_high, sizeof(double)); + V_rho = 4/(a*a*a); V_my_s = (V_rho * 100 * sigma_inc); V_my_a_v = (V_rho * 100 * sigma_abs * 2200); @@ -403,6 +401,8 @@ INITIALIZE %} TRACE %{ + double* vf_list = (double *)calloc( e_steps_low*e_steps_high, sizeof(double)); // List of allowed final velocities. Has length of scan_steps + struct neutron_params neutron; double t0, t1; /* Entry/exit time for cylinder */ double v_i, v_f; /* Neutron velocities: initial, final */ double vx_i, vy_i, vz_i; /* Neutron initial velocity vector */ @@ -502,16 +502,14 @@ TRACE p4 = 2*VS2E*v_f/J_factor; /* Jacobian of delta functions in cross section */ p5 = b*b/M; /* Cross section factor 2 */ p *= p1*p2*p3*p4*p5; - // Finally, set vf_list to zero again - memset(vf_list, 0, nf * sizeof(double)); } /* else transmit: Neutron did not hit the sample */ - +free(vf_list); %} FINALLY %{ - free(vf_list); + %} From 41c541ea4f9fa5957997c0f4352aba9d5d8e957a Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 20 Aug 2025 14:12:29 +0200 Subject: [PATCH 6/7] Primary: Fixes to findroot_gpu function (sync from base findroot) to allow correct GPU execution. Secondary changes: 1. calloc not available on gpu (in TRACE), use malloc 2. Don't use e_steps_low*e_steps_high but e_steps_low+e_steps_high for array allocation 3. Cosmetics to comp header --- mcstas-comps/samples/Phonon_simple.comp | 53 ++++++++++++++++++------- 1 file changed, 38 insertions(+), 15 deletions(-) diff --git a/mcstas-comps/samples/Phonon_simple.comp b/mcstas-comps/samples/Phonon_simple.comp index e84fb1f2a7..0675a116f2 100644 --- a/mcstas-comps/samples/Phonon_simple.comp +++ b/mcstas-comps/samples/Phonon_simple.comp @@ -51,8 +51,8 @@ * target_z: [m] position of target to focus at. Straight ahead. * target_index: [1] relative index of component to focus at, e.g. next is +1 * gap: [meV] Bandgap energy (unphysical) -* e_steps_low [50] Amount of possible intersections beneath the elastic line -* e_steps_high [50] Amount of possible intersections above the elastic line +* e_steps_low [1] Amount of possible intersections beneath the elastic line +* e_steps_high [1] Amount of possible intersections above the elastic line * * * CALCULATED PARAMETERS: @@ -307,7 +307,7 @@ double zridd_gpu(double x1, double x2, struct neutron_params *neutron, struct ph double root; // Energy gain and energy loss spaces are not equally big. We check uniformly // So we use two different ranges - double range_low = brack_mid - brack_low; + double range_low = brack_mid - brack_low; double range_high = brack_high - brack_mid; // First in energy loss for the neutron for (int i=0; ie_steps_low_; i++){ @@ -335,25 +335,33 @@ double zridd_gpu(double x1, double x2, struct neutron_params *neutron, struct ph void findroots_gpu(double brack_low, double brack_mid, double brack_high, double *list, int* index, struct neutron_params *neutron, struct phonon_params *phonon) { - double root,range=brack_mid-brack_low; - int i, steps=100; - - for (i=0; ie_steps_low_; i++){ + root = zridd_gpu(brack_low+range_low*i/ phonon->e_steps_low_, + brack_low+range_low*(i+1)/ phonon->e_steps_low_, + neutron, phonon, ROOTACC); if (root != UNUSED) { list[(*index)++]=root; } } - root = zridd_gpu(brack_mid, brack_high, neutron,phonon, ROOTACC); - if (root != UNUSED) - { + + // Then in energy gain for the neutron + for (int i=0; ie_steps_high_; i++){ + root = zridd_gpu(brack_low+range_high*i/ phonon->e_steps_high_, + brack_low+range_high*(i+1)/ phonon->e_steps_high_, + neutron, phonon, ROOTACC); + if (root != UNUSED){ list[(*index)++]=root; } + } } + #undef UNUSED #undef MAXRIDD @@ -401,7 +409,22 @@ INITIALIZE %} TRACE %{ - double* vf_list = (double *)calloc( e_steps_low*e_steps_high, sizeof(double)); // List of allowed final velocities. Has length of scan_steps + double* vf_list; + #ifdef OPENACC + vf_list = (double *)malloc( (e_steps_low + e_steps_high)* sizeof(double)); // List of allowed final velocities. Has length of scan_steps + #else + vf_list = (double *)calloc( e_steps_low + e_steps_high, sizeof(double)); // List of allowed final velocities. Has length of scan_steps + #endif + if (!vf_list) { + printf("Memory allocation failed, fatal error!\n"); + exit(-1); + } + #ifdef OPENACC + for (int ii=0; ii Date: Wed, 20 Aug 2025 14:21:27 +0200 Subject: [PATCH 7/7] Purely indentation fixes / formatting --- mcstas-comps/samples/Phonon_simple.comp | 69 ++++++++++++------------- 1 file changed, 32 insertions(+), 37 deletions(-) diff --git a/mcstas-comps/samples/Phonon_simple.comp b/mcstas-comps/samples/Phonon_simple.comp index 0675a116f2..29a1e88750 100644 --- a/mcstas-comps/samples/Phonon_simple.comp +++ b/mcstas-comps/samples/Phonon_simple.comp @@ -302,67 +302,62 @@ double zridd_gpu(double x1, double x2, struct neutron_params *neutron, struct ph void findroots(double brack_low, double brack_mid, double brack_high, double *list, int* index, double (*f)(struct neutron_params*, struct phonon_params*), - struct neutron_params *neutron, struct phonon_params *phonon) - { - double root; - // Energy gain and energy loss spaces are not equally big. We check uniformly - // So we use two different ranges - double range_low = brack_mid - brack_low; - double range_high = brack_high - brack_mid; - // First in energy loss for the neutron - for (int i=0; ie_steps_low_; i++){ + struct neutron_params *neutron, struct phonon_params *phonon) { + double root; + // Energy gain and energy loss spaces are not equally big. We check uniformly + // So we use two different ranges + double range_low = brack_mid - brack_low; + double range_high = brack_high - brack_mid; + // First in energy loss for the neutron + for (int i=0; ie_steps_low_; i++){ root = zridd(f, brack_low+range_low*i/ phonon->e_steps_low_, brack_low+range_low*(i+1)/ phonon->e_steps_low_, neutron, phonon, ROOTACC); - if (root != UNUSED) - { - list[(*index)++]=root; + if (root != UNUSED) { + list[(*index)++]=root; } - } + } // Then in energy gain for the neutron for (int i=0; ie_steps_high_; i++){ root = zridd(f, brack_low+range_high*i/ phonon->e_steps_high_, - brack_low+range_high*(i+1)/ phonon->e_steps_high_, - neutron, phonon, ROOTACC); + brack_low+range_high*(i+1)/ phonon->e_steps_high_, + neutron, phonon, ROOTACC); if (root != UNUSED){ - list[(*index)++]=root; + list[(*index)++]=root; } - } } + } #pragma acc routine void findroots_gpu(double brack_low, double brack_mid, double brack_high, - double *list, int* index, struct neutron_params *neutron, struct phonon_params *phonon) - { - double root; - // Energy gain and energy loss spaces are not equally big. We check uniformly - // So we use two different ranges - double range_low = brack_mid - brack_low; - double range_high = brack_high - brack_mid; - // First in energy loss for the neutron - for (int i=0; ie_steps_low_; i++){ + double *list, int* index, struct neutron_params *neutron, struct phonon_params *phonon) { + double root; + // Energy gain and energy loss spaces are not equally big. We check uniformly + // So we use two different ranges + double range_low = brack_mid - brack_low; + double range_high = brack_high - brack_mid; + // First in energy loss for the neutron + for (int i=0; ie_steps_low_; i++){ root = zridd_gpu(brack_low+range_low*i/ phonon->e_steps_low_, - brack_low+range_low*(i+1)/ phonon->e_steps_low_, - neutron, phonon, ROOTACC); - if (root != UNUSED) - { - list[(*index)++]=root; + brack_low+range_low*(i+1)/ phonon->e_steps_low_, + neutron, phonon, ROOTACC); + if (root != UNUSED) { + list[(*index)++]=root; } - } + } // Then in energy gain for the neutron for (int i=0; ie_steps_high_; i++){ root = zridd_gpu(brack_low+range_high*i/ phonon->e_steps_high_, - brack_low+range_high*(i+1)/ phonon->e_steps_high_, - neutron, phonon, ROOTACC); + brack_low+range_high*(i+1)/ phonon->e_steps_high_, + neutron, phonon, ROOTACC); if (root != UNUSED){ - list[(*index)++]=root; + list[(*index)++]=root; } - } } + } - #undef UNUSED #undef MAXRIDD #endif