Changeset 2c108a3 in sasmodels


Ignore:
Timestamp:
Oct 19, 2017 12:58:57 PM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
6773b02
Parents:
9ee2756
Message:

simplify magnetism code

Files:
4 edited

Legend:

Unmodified
Added
Removed
  • doc/developer/calculator.rst

    r9ee2756 r2c108a3  
    107107- MAGNETIC_PARS is a comma separated list of the magnetic SLDs, indicating 
    108108  their locations in the values vector. 
    109 - MAGNETIC_PAR1, ... are the first three magnetic parameter ids so we can 
    110   hard code the setting of magnetic values if there are only a few of them. 
    111109- KERNEL_NAME is the name of the function being declared 
    112110- PARAMETER_TABLE is the declaration of the parameters to the kernel: 
     
    208206    NUM_MAGNETIC = 2      // two parameters might be magnetic 
    209207    MAGNETIC_PARS = 4, 5  // they are sld and sld_solvent 
    210     MAGNETIC_PAR0 = 4     // sld index 
    211     MAGNETIC_PAR1 = 5     // solvent index 
    212208 
    213209    details { 
  • doc/guide/magnetism/magnetism.rst

    r1f058ea r2c108a3  
    3131to the $x'$ axis, the possible spin states after the sample are then 
    3232 
    33 No spin-flips (+ +) and (- -) 
     33Non spin-flip (+ +) and (- -) 
    3434 
    35 Spin-flips    (+ -) and (- +) 
     35Spin-flip    (+ -) and (- +) 
    3636 
    3737.. figure:: 
     
    4141$\phi$ and $\theta_{up}$, respectively, then, depending on the spin state of the 
    4242neutrons, the scattering length densities, including the nuclear scattering 
    43 length density ($\beta{_N}$) are 
     43length density $(\beta{_N})$ are 
    4444 
    4545.. math:: 
    4646    \beta_{\pm\pm} =  \beta_N \mp D_M M_{\perp x'} 
    47     \text{ when there are no spin-flips} 
     47    \text{ for non spin-flip states} 
    4848 
    4949and 
     
    5151.. math:: 
    5252    \beta_{\pm\mp} =  -D_M (M_{\perp y'} \pm iM_{\perp z'}) 
    53     \text{ when there are} 
     53    \text{ for spin-flip states} 
    5454 
    5555where 
    5656 
    5757.. math:: 
    58     M_{\perp x'} = M_{0q_x}\cos(\theta_{up})+M_{0q_y}\sin(\theta_{up}) \\ 
    59     M_{\perp y'} = M_{0q_y}\cos(\theta_{up})-M_{0q_x}\sin(\theta_{up}) \\ 
    60     M_{\perp z'} = M_{0z} \\ 
    61     M_{0q_x} = (M_{0x}\cos\phi - M_{0y}\sin\phi)\cos\phi \\ 
    62     M_{0q_y} = (M_{0y}\sin\phi - M_{0x}\cos\phi)\sin\phi 
     58    M_{\perp x'} &= M_{0q_x}\cos(\theta_{up})+M_{0q_y}\sin(\theta_{up}) \\ 
     59    M_{\perp y'} &= M_{0q_y}\cos(\theta_{up})-M_{0q_x}\sin(\theta_{up}) \\ 
     60    M_{\perp z'} &= M_{0z} \\ 
     61    M_{0q_x} &= (M_{0x}\cos\phi - M_{0y}\sin\phi)\cos\phi \\ 
     62    M_{0q_y} &= (M_{0y}\sin\phi - M_{0x}\cos\phi)\sin\phi 
    6363 
    6464Here, $M_{0x}$, $M_{0x}$, $M_{0z}$ are the x, y and z components 
     
    6666 
    6767.. math:: 
    68     M_{0x} = M_0\cos\theta_M\cos\phi_M \\ 
    69     M_{0y} = M_0\sin\theta_M \\ 
    70     M_{0z} = -M_0\cos\theta_M\sin\phi_M 
     68    M_{0x} &= M_0\cos\theta_M\cos\phi_M \\ 
     69    M_{0y} &= M_0\sin\theta_M \\ 
     70    M_{0z} &= -M_0\cos\theta_M\sin\phi_M 
    7171 
    7272and the magnetization angles $\theta_M$ and $\phi_M$ are defined in 
  • sasmodels/generate.py

    r9ee2756 r2c108a3  
    760760    source.append("#define NUM_MAGNETIC %d" % partable.nmagnetic) 
    761761    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    762     for k, v in enumerate(magpars[:3]): 
    763         source.append("#define MAGNETIC_PAR%d %d"%(k+1, v)) 
    764762 
    765763    # TODO: allow mixed python/opencl kernels? 
  • sasmodels/kernel_iq.c

    r9ee2756 r2c108a3  
    1717//  NUM_VALUES : the number of values to skip at the start of the 
    1818//      values array before you get to the dispersity values. 
    19 //  NUM_MAGNETIC : the number of magnetic parameters 
    2019//  PARAMETER_TABLE : list of parameter declarations used to create the 
    2120//      ParameterTable type. 
    2221//  KERNEL_NAME : model_Iq, model_Iqxy or model_Imagnetic.  This code is 
    2322//      included three times, once for each kernel type. 
     23//  MAGNETIC : defined when the magnetic kernel is being instantiated 
     24//  NUM_MAGNETIC : the number of magnetic parameters 
    2425//  MAGNETIC_PARS : a comma-separated list of indices to the sld 
    2526//      parameters in the parameter table. 
    26 //  MAGNETIC_PAR1, ... : the first few indices of slds in the parameter table. 
    27 //  MAGNETIC : defined when the magnetic kernel is being instantiated 
    2827//  CALL_VOLUME(table) : call the form volume function 
    2928//  CALL_IQ(q, table) : call the Iq function for 1D calcs. 
     
    3938typedef struct { 
    4039#if MAX_PD > 0 
    41     int32_t pd_par[MAX_PD];     // id of the nth polydispersity variable 
    42     int32_t pd_length[MAX_PD];  // length of the nth polydispersity weight vector 
     40    int32_t pd_par[MAX_PD];     // id of the nth dispersity variable 
     41    int32_t pd_length[MAX_PD];  // length of the nth dispersity weight vector 
    4342    int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector 
    4443    int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level 
     
    7372//     uu * (sld - m_sigma_x); 
    7473//     dd * (sld + m_sigma_x); 
    75 //     ud * (m_sigma_y + 1j*m_sigma_z); 
    76 //     du * (m_sigma_y - 1j*m_sigma_z); 
    77 static void set_spins(double in_spin, double out_spin, double spins[4]) 
     74//     ud * (m_sigma_y - 1j*m_sigma_z); 
     75//     du * (m_sigma_y + 1j*m_sigma_z); 
     76// weights for spin crosssections: dd du real, ud real, uu, du imag, ud imag 
     77static void set_spin_weights(double in_spin, double out_spin, double spins[4]) 
    7878{ 
    7979  in_spin = clip(in_spin, 0.0, 1.0); 
    8080  out_spin = clip(out_spin, 0.0, 1.0); 
    8181  spins[0] = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); // dd 
    82   spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du 
    83   spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud 
     82  spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du real 
     83  spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud real 
    8484  spins[3] = sqrt(sqrt(in_spin * out_spin));             // uu 
     85  spins[4] = spins[1]; // du imag 
     86  spins[5] = spins[2]; // ud imag 
    8587} 
    8688 
    8789// Compute the magnetic sld 
    88 static double mag_sld(double qx, double qy, double p, 
    89                        double mx, double my, double sld) 
    90 { 
     90static double mag_sld( 
     91  const int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag 
     92  const double qx, const double qy, 
     93  const double px, const double py, 
     94  const double sld, 
     95  const double mx, const double my, const double mz 
     96) 
     97{ 
     98  if (xs < 4) { 
    9199    const double perp = qy*mx - qx*my; 
    92     return sld + perp*p; 
    93 } 
     100    switch (xs) { 
     101      case 0: // uu => sld - D M_perpx 
     102          return sld - px*perp; 
     103      case 1: // ud real => -D M_perpy 
     104          return py*perp; 
     105      case 2: // du real => -D M_perpy 
     106          return py*perp; 
     107      case 3: // dd real => sld + D M_perpx 
     108          return sld + px*perp; 
     109    } 
     110  } else { 
     111    if (xs== 4) { 
     112      return -mz;  // ud imag => -D M_perpz 
     113    } else { // index == 5 
     114      return mz;   // du imag => D M_perpz 
     115    } 
     116  } 
     117} 
     118 
    94119 
    95120#endif 
     
    235260void KERNEL_NAME( 
    236261    int32_t nq,                 // number of q values 
    237     const int32_t pd_start,     // where we are in the polydispersity loop 
    238     const int32_t pd_stop,      // where we are stopping in the polydispersity loop 
     262    const int32_t pd_start,     // where we are in the dispersity loop 
     263    const int32_t pd_stop,      // where we are stopping in the dispersity loop 
    239264    global const ProblemDetails *details, 
    240265    global const double *values, 
    241266    global const double *q, // nq q values, with padding to boundary 
    242267    global double *result,  // nq+1 return values, again with padding 
    243     const double cutoff     // cutoff in the polydispersity weight product 
     268    const double cutoff     // cutoff in the dispersity weight product 
    244269    ) 
    245270{ 
     
    256281 
    257282  // Storage for the current parameter values.  These will be updated as we 
    258   // walk the polydispersity cube. 
     283  // walk the dispersity mesh. 
    259284  ParameterBlock local_values; 
    260285 
     
    262287  // Location of the sld parameters in the parameter vector. 
    263288  // These parameters are updated with the effective sld due to magnetism. 
    264   #if NUM_MAGNETIC > 3 
    265289  const int32_t slds[] = { MAGNETIC_PARS }; 
    266   #endif 
    267290 
    268291  // Interpret polarization cross section. 
     
    271294  //     up_angle = values[NUM_PARS+4]; 
    272295  // TODO: could precompute more magnetism parameters before calling the kernel. 
    273   double spins[4]; 
     296  double spins[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill 
    274297  double cos_mspin, sin_mspin; 
    275   set_spins(values[NUM_PARS+2], values[NUM_PARS+3], spins); 
     298  set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], spins); 
    276299  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin); 
    277300#endif // MAGNETIC 
     
    395418  int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP; 
    396419 
    397 // Jump into the middle of the polydispersity loop 
     420// Jump into the middle of the dispersity loop 
    398421#define PD_OPEN(_LOOP,_OUTER) \ 
    399422  while (i##_LOOP < n##_LOOP) { \ 
     
    427450  // unoriented 1D 
    428451  double qk; 
    429   #define FETCH_Q do { qk = q[q_index]; } while (0) 
    430   #define BUILD_ROTATION do {} while(0) 
    431   #define APPLY_ROTATION do {} while(0) 
    432   #define CALL_KERNEL CALL_IQ(qk, local_values.table) 
     452  #define FETCH_Q() do { qk = q[q_index]; } while (0) 
     453  #define BUILD_ROTATION() do {} while(0) 
     454  #define APPLY_ROTATION() do {} while(0) 
     455  #define CALL_KERNEL() CALL_IQ(qk, local_values.table) 
    433456 
    434457#elif defined(CALL_IQ_A) 
    435458  // unoriented 2D 
    436459  double qx, qy; 
    437   #define FETCH_Q do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
    438   #define BUILD_ROTATION do {} while(0) 
    439   #define APPLY_ROTATION do {} while(0) 
    440   #define CALL_KERNEL CALL_IQ_A(sqrt(qx*qx+qy*qy), local_values.table) 
     460  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
     461  #define BUILD_ROTATION() do {} while(0) 
     462  #define APPLY_ROTATION() do {} while(0) 
     463  #define CALL_KERNEL() CALL_IQ_A(sqrt(qx*qx+qy*qy), local_values.table) 
    441464 
    442465#elif defined(CALL_IQ_AC) 
    443466  // oriented symmetric 2D 
    444467  double qx, qy; 
    445   #define FETCH_Q do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
     468  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
    446469  double qa, qc; 
    447470  QACRotation rotation; 
     
    452475  const double theta = values[details->theta_par+2]; 
    453476  const double phi = values[details->theta_par+3]; 
    454   #define BUILD_ROTATION qac_rotation(&rotation, \ 
     477  #define BUILD_ROTATION() qac_rotation(&rotation, \ 
    455478    theta, phi, local_values.table.theta, local_values.table.phi) 
    456   #define APPLY_ROTATION qac_apply(rotation, qx, qy, &qa, &qc) 
    457   #define CALL_KERNEL CALL_IQ_AC(qa, qc, local_values.table) 
     479  #define APPLY_ROTATION() qac_apply(rotation, qx, qy, &qa, &qc) 
     480  #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table) 
    458481 
    459482#elif defined(CALL_IQ_ABC) 
    460483  // oriented asymmetric 2D 
    461484  double qx, qy; 
    462   #define FETCH_Q do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
     485  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
    463486  double qa, qb, qc; 
    464487  QABCRotation rotation; 
     
    470493  const double phi = values[details->theta_par+3]; 
    471494  const double psi = values[details->theta_par+4]; 
    472   #define BUILD_ROTATION qabc_rotation(&rotation, \ 
     495  #define BUILD_ROTATION() qabc_rotation(&rotation, \ 
    473496    theta, phi, psi, local_values.table.theta, \ 
    474497    local_values.table.phi, local_values.table.psi) 
    475   #define APPLY_ROTATION qabc_apply(rotation, qx, qy, &qa, &qb, &qc) 
    476   #define CALL_KERNEL CALL_IQ_ABC(qa, qb, qc, local_values.table) 
     498  #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc) 
     499  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 
    477500 
    478501#endif 
     
    538561    if (weight0 > cutoff) { 
    539562      pd_norm += weight0 * CALL_VOLUME(local_values.table); 
    540       BUILD_ROTATION; 
     563      BUILD_ROTATION(); 
    541564 
    542565#ifndef USE_OPENCL 
     
    549572      { 
    550573 
    551         FETCH_Q; 
    552         APPLY_ROTATION; 
     574        FETCH_Q(); 
     575        APPLY_ROTATION(); 
    553576 
    554577        // ======= COMPUTE SCATTERING ========== 
    555578        #if defined(MAGNETIC) && NUM_MAGNETIC > 0 
    556         // Constant across orientation, polydispersity for given qx, qy 
    557         double scattering = 0.0; 
    558         // TODO: what is the magnetic scattering at q=0 
    559         const double qsq = qx*qx + qy*qy; 
    560         if (qsq > 1.e-16) { 
    561           double p[4];  // dd, du, ud, uu 
    562           p[0] = (qy*cos_mspin + qx*sin_mspin)/qsq; 
    563           p[3] = -p[0]; 
    564           p[1] = p[2] = (qy*sin_mspin - qx*cos_mspin)/qsq; 
    565  
    566           for (int index=0; index<4; index++) { 
    567             const double xs = spins[index]; 
    568             if (xs > 1.e-8) { 
    569               const int spin_flip = (index==1) || (index==2); 
    570               const double pk = p[index]; 
    571               for (int axis=0; axis<=spin_flip; axis++) { 
    572                 #define SLD(_M_offset, _sld_offset) \ 
    573                     local_values.vector[_sld_offset] = xs * (axis \ 
    574                     ? (index==1 ? -values[_M_offset+2] : values[_M_offset+2]) \ 
    575                     : mag_sld(qx, qy, pk, values[_M_offset], values[_M_offset+1], \ 
    576                               (spin_flip ? 0.0 : values[_sld_offset+2]))) 
    577                 #define M1 NUM_PARS+5 
    578                 #if NUM_MAGNETIC==1 
    579                   SLD(M1, MAGNETIC_PAR1); 
    580                 #elif NUM_MAGNETIC==2 
    581                   SLD(M1, MAGNETIC_PAR1); 
    582                   SLD(M1+3, MAGNETIC_PAR2); 
    583                 #elif NUM_MAGNETIC==3 
    584                   SLD(M1, MAGNETIC_PAR1); 
    585                   SLD(M1+3, MAGNETIC_PAR2); 
    586                   SLD(M1+6, MAGNETIC_PAR3); 
    587                 #else 
    588                   for (int sk=0; sk<NUM_MAGNETIC; sk++) { 
    589                       SLD(M1+3*sk, slds[sk]); 
    590                   } 
    591                 #endif 
    592                 #undef SLD 
    593                 scattering += CALL_KERNEL; 
     579          // Compute the scattering from the magnetic cross sections. 
     580          double scattering = 0.0; 
     581          const double qsq = qx*qx + qy*qy; 
     582          if (qsq > 1.e-16) { 
     583            // TODO: what is the magnetic scattering at q=0 
     584            const double px = (qy*cos_mspin + qx*sin_mspin)/qsq; 
     585            const double py = (qy*sin_mspin - qx*cos_mspin)/qsq; 
     586 
     587            // loop over uu, ud real, du real, dd, ud imag, du imag 
     588            for (int xs=0; xs<6; xs++) { 
     589              const double xs_weight = spins[xs]; 
     590              if (xs_weight > 1.e-8) { 
     591                // Since the cross section weight is significant, set the slds 
     592                // to the effective slds for this cross section, call the 
     593                // kernel, and add according to weight. 
     594                for (int sk=0; sk<NUM_MAGNETIC; sk++) { 
     595                  const int32_t mag_index = NUM_PARS+5 + 3*sk; 
     596                  const int32_t sld_index = slds[sk]; 
     597                  const double mx = values[mag_index]; 
     598                  const double my = values[mag_index+1]; 
     599                  const double mz = values[mag_index+2]; 
     600                  local_values.vector[sld_index] = 
     601                    mag_sld(xs, qx, qy, px, py, values[sld_index+2], mx, my, mz); 
     602                } 
     603                scattering += xs_weight * CALL_KERNEL(); 
    594604              } 
    595605            } 
    596606          } 
    597         } 
    598607        #else  // !MAGNETIC 
    599         const double scattering = CALL_KERNEL; 
     608          const double scattering = CALL_KERNEL(); 
    600609        #endif // !MAGNETIC 
    601610//printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight, weight0); 
    602611 
    603612        #ifdef USE_OPENCL 
    604         this_result += weight0 * scattering; 
     613          this_result += weight0 * scattering; 
    605614        #else // !USE_OPENCL 
    606         result[q_index] += weight0 * scattering; 
     615          result[q_index] += weight0 * scattering; 
    607616        #endif // !USE_OPENCL 
    608617      } 
     
    627636  PD_CLOSE(4) 
    628637#endif 
    629  
    630 // clear the macros in preparation for the next kernel 
    631 #undef PD_INIT 
    632 #undef PD_OPEN 
    633 #undef PD_CLOSE 
    634 #undef FETCH_Q 
    635 #undef BUILD_ROTATION 
    636 #undef APPLY_ROTATION 
    637 #undef CALL_KERNEL 
    638638 
    639639// Remember the current result and the updated norm. 
     
    646646//printf("res: %g/%g\n", result[0], pd_norm); 
    647647#endif // !USE_OPENCL 
    648 } 
     648 
     649// clear the macros in preparation for the next kernel 
     650#undef PD_INIT 
     651#undef PD_OPEN 
     652#undef PD_CLOSE 
     653#undef FETCH_Q 
     654#undef BUILD_ROTATION 
     655#undef APPLY_ROTATION 
     656#undef CALL_KERNEL 
     657} 
Note: See TracChangeset for help on using the changeset viewer.