Changeset 2c108a3 in sasmodels

Timestamp:

Oct 19, 2017 2:58:57 PM (7 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:

• doc/developer/calculator.rst (modified) (2 diffs)
• doc/guide/magnetism/magnetism.rst (modified) (4 diffs)
• sasmodels/generate.py (modified) (1 diff)
• sasmodels/kernel_iq.c (modified) (15 diffs)

doc/developer/calculator.rst

@@ -107 +107 @@
 - MAGNETIC_PARS is a comma separated list of the magnetic SLDs, indicating
   their locations in the values vector.
-- MAGNETIC_PAR1, ... are the first three magnetic parameter ids so we can
-  hard code the setting of magnetic values if there are only a few of them.
 - KERNEL_NAME is the name of the function being declared
 - PARAMETER_TABLE is the declaration of the parameters to the kernel:
@@ -208 +206 @@
     NUM_MAGNETIC = 2      // two parameters might be magnetic
     MAGNETIC_PARS = 4, 5  // they are sld and sld_solvent
-    MAGNETIC_PAR0 = 4     // sld index
-    MAGNETIC_PAR1 = 5     // solvent index
 
     details {

doc/guide/magnetism/magnetism.rst

@@ -31 +31 @@
 to the $x'$ axis, the possible spin states after the sample are then
 
-No spin-flips (+ +) and (- -)
+Non spin-flip (+ +) and (- -)
 
-Spin-flips    (+ -) and (- +)
+Spin-flip    (+ -) and (- +)
 
 .. figure::
@@ -41 +41 @@
 $\phi$ and $\theta_{up}$, respectively, then, depending on the spin state of the
 neutrons, the scattering length densities, including the nuclear scattering
-length density ($\beta{_N}$) are
+length density $(\beta{_N})$ are
 
 .. math::
     \beta_{\pm\pm} =  \beta_N \mp D_M M_{\perp x'}
-    \text{ when there are no spin-flips}
+    \text{ for non spin-flip states}
 
 and
@@ -51 +51 @@
 .. math::
     \beta_{\pm\mp} =  -D_M (M_{\perp y'} \pm iM_{\perp z'})
-    \text{ when there are}
+    \text{ for spin-flip states}
 
 where
 
 .. math::
-    M_{\perp x'} = M_{0q_x}\cos(\theta_{up})+M_{0q_y}\sin(\theta_{up}) \\
-    M_{\perp y'} = M_{0q_y}\cos(\theta_{up})-M_{0q_x}\sin(\theta_{up}) \\
-    M_{\perp z'} = M_{0z} \\
-    M_{0q_x} = (M_{0x}\cos\phi - M_{0y}\sin\phi)\cos\phi \\
-    M_{0q_y} = (M_{0y}\sin\phi - M_{0x}\cos\phi)\sin\phi
+    M_{\perp x'} &= M_{0q_x}\cos(\theta_{up})+M_{0q_y}\sin(\theta_{up}) \\
+    M_{\perp y'} &= M_{0q_y}\cos(\theta_{up})-M_{0q_x}\sin(\theta_{up}) \\
+    M_{\perp z'} &= M_{0z} \\
+    M_{0q_x} &= (M_{0x}\cos\phi - M_{0y}\sin\phi)\cos\phi \\
+    M_{0q_y} &= (M_{0y}\sin\phi - M_{0x}\cos\phi)\sin\phi
 
 Here, $M_{0x}$, $M_{0x}$, $M_{0z}$ are the x, y and z components
@@ -66 +66 @@
 
 .. math::
-    M_{0x} = M_0\cos\theta_M\cos\phi_M \\
-    M_{0y} = M_0\sin\theta_M \\
-    M_{0z} = -M_0\cos\theta_M\sin\phi_M
+    M_{0x} &= M_0\cos\theta_M\cos\phi_M \\
+    M_{0y} &= M_0\sin\theta_M \\
+    M_{0z} &= -M_0\cos\theta_M\sin\phi_M
 
 and the magnetization angles $\theta_M$ and $\phi_M$ are defined in

sasmodels/generate.py

@@ -760 +760 @@
     source.append("#define NUM_MAGNETIC %d" % partable.nmagnetic)
     source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars))
-    for k, v in enumerate(magpars[:3]):
-        source.append("#define MAGNETIC_PAR%d %d"%(k+1, v))
 
     # TODO: allow mixed python/opencl kernels?

sasmodels/kernel_iq.c

@@ -17 +17 @@
 //  NUM_VALUES : the number of values to skip at the start of the
 //      values array before you get to the dispersity values.
-//  NUM_MAGNETIC : the number of magnetic parameters
 //  PARAMETER_TABLE : list of parameter declarations used to create the
 //      ParameterTable type.
 //  KERNEL_NAME : model_Iq, model_Iqxy or model_Imagnetic.  This code is
 //      included three times, once for each kernel type.
+//  MAGNETIC : defined when the magnetic kernel is being instantiated
+//  NUM_MAGNETIC : the number of magnetic parameters
 //  MAGNETIC_PARS : a comma-separated list of indices to the sld
 //      parameters in the parameter table.
-//  MAGNETIC_PAR1, ... : the first few indices of slds in the parameter table.
-//  MAGNETIC : defined when the magnetic kernel is being instantiated
 //  CALL_VOLUME(table) : call the form volume function
 //  CALL_IQ(q, table) : call the Iq function for 1D calcs.
@@ -39 +38 @@
 typedef struct {
 #if MAX_PD > 0
-    int32_t pd_par[MAX_PD];     // id of the nth polydispersity variable
-    int32_t pd_length[MAX_PD];  // length of the nth polydispersity weight vector
+    int32_t pd_par[MAX_PD];     // id of the nth dispersity variable
+    int32_t pd_length[MAX_PD];  // length of the nth dispersity weight vector
     int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector
     int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level
@@ -73 +72 @@
 //     uu * (sld - m_sigma_x);
 //     dd * (sld + m_sigma_x);
-//     ud * (m_sigma_y + 1j*m_sigma_z);
-//     du * (m_sigma_y - 1j*m_sigma_z);
-static void set_spins(double in_spin, double out_spin, double spins[4])
+//     ud * (m_sigma_y - 1j*m_sigma_z);
+//     du * (m_sigma_y + 1j*m_sigma_z);
+// weights for spin crosssections: dd du real, ud real, uu, du imag, ud imag
+static void set_spin_weights(double in_spin, double out_spin, double spins[4])
 {
   in_spin = clip(in_spin, 0.0, 1.0);
   out_spin = clip(out_spin, 0.0, 1.0);
   spins[0] = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); // dd
-  spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du
-  spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud
+  spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du real
+  spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud real
   spins[3] = sqrt(sqrt(in_spin * out_spin));             // uu
+  spins[4] = spins[1]; // du imag
+  spins[5] = spins[2]; // ud imag
 }
 
 // Compute the magnetic sld
-static double mag_sld(double qx, double qy, double p,
-                       double mx, double my, double sld)
-{
+static double mag_sld(
+  const int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag
+  const double qx, const double qy,
+  const double px, const double py,
+  const double sld,
+  const double mx, const double my, const double mz
+)
+{
+  if (xs < 4) {
     const double perp = qy*mx - qx*my;
-    return sld + perp*p;
-}
+    switch (xs) {
+      case 0: // uu => sld - D M_perpx
+          return sld - px*perp;
+      case 1: // ud real => -D M_perpy
+          return py*perp;
+      case 2: // du real => -D M_perpy
+          return py*perp;
+      case 3: // dd real => sld + D M_perpx
+          return sld + px*perp;
+    }
+  } else {
+    if (xs== 4) {
+      return -mz;  // ud imag => -D M_perpz
+    } else { // index == 5
+      return mz;   // du imag => D M_perpz
+    }
+  }
+}
+
 
 #endif
@@ -235 +260 @@
 void KERNEL_NAME(
     int32_t nq,                 // number of q values
-    const int32_t pd_start,     // where we are in the polydispersity loop
-    const int32_t pd_stop,      // where we are stopping in the polydispersity loop
+    const int32_t pd_start,     // where we are in the dispersity loop
+    const int32_t pd_stop,      // where we are stopping in the dispersity loop
     global const ProblemDetails *details,
     global const double *values,
     global const double *q, // nq q values, with padding to boundary
     global double *result,  // nq+1 return values, again with padding
-    const double cutoff     // cutoff in the polydispersity weight product
+    const double cutoff     // cutoff in the dispersity weight product
     )
 {
@@ -256 +281 @@
 
   // Storage for the current parameter values.  These will be updated as we
-  // walk the polydispersity cube.
+  // walk the dispersity mesh.
   ParameterBlock local_values;
 
@@ -262 +287 @@
   // Location of the sld parameters in the parameter vector.
   // These parameters are updated with the effective sld due to magnetism.
-  #if NUM_MAGNETIC > 3
   const int32_t slds[] = { MAGNETIC_PARS };
-  #endif
 
   // Interpret polarization cross section.
@@ -271 +294 @@
   //     up_angle = values[NUM_PARS+4];
   // TODO: could precompute more magnetism parameters before calling the kernel.
-  double spins[4];
+  double spins[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill
   double cos_mspin, sin_mspin;
-  set_spins(values[NUM_PARS+2], values[NUM_PARS+3], spins);
+  set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], spins);
   SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin);
 #endif // MAGNETIC
@@ -395 +418 @@
   int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP;
 
-// Jump into the middle of the polydispersity loop
+// Jump into the middle of the dispersity loop
 #define PD_OPEN(_LOOP,_OUTER) \
   while (i##_LOOP < n##_LOOP) { \
@@ -427 +450 @@
   // unoriented 1D
   double qk;
-  #define FETCH_Q do { qk = q[q_index]; } while (0)
-  #define BUILD_ROTATION do {} while(0)
-  #define APPLY_ROTATION do {} while(0)
-  #define CALL_KERNEL CALL_IQ(qk, local_values.table)
+  #define FETCH_Q() do { qk = q[q_index]; } while (0)
+  #define BUILD_ROTATION() do {} while(0)
+  #define APPLY_ROTATION() do {} while(0)
+  #define CALL_KERNEL() CALL_IQ(qk, local_values.table)
 
 #elif defined(CALL_IQ_A)
   // unoriented 2D
   double qx, qy;
-  #define FETCH_Q do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
-  #define BUILD_ROTATION do {} while(0)
-  #define APPLY_ROTATION do {} while(0)
-  #define CALL_KERNEL CALL_IQ_A(sqrt(qx*qx+qy*qy), local_values.table)
+  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
+  #define BUILD_ROTATION() do {} while(0)
+  #define APPLY_ROTATION() do {} while(0)
+  #define CALL_KERNEL() CALL_IQ_A(sqrt(qx*qx+qy*qy), local_values.table)
 
 #elif defined(CALL_IQ_AC)
   // oriented symmetric 2D
   double qx, qy;
-  #define FETCH_Q do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
+  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
   double qa, qc;
   QACRotation rotation;
@@ -452 +475 @@
   const double theta = values[details->theta_par+2];
   const double phi = values[details->theta_par+3];
-  #define BUILD_ROTATION qac_rotation(&rotation, \
+  #define BUILD_ROTATION() qac_rotation(&rotation, \
     theta, phi, local_values.table.theta, local_values.table.phi)
-  #define APPLY_ROTATION qac_apply(rotation, qx, qy, &qa, &qc)
-  #define CALL_KERNEL CALL_IQ_AC(qa, qc, local_values.table)
+  #define APPLY_ROTATION() qac_apply(rotation, qx, qy, &qa, &qc)
+  #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table)
 
 #elif defined(CALL_IQ_ABC)
   // oriented asymmetric 2D
   double qx, qy;
-  #define FETCH_Q do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
+  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
   double qa, qb, qc;
   QABCRotation rotation;
@@ -470 +493 @@
   const double phi = values[details->theta_par+3];
   const double psi = values[details->theta_par+4];
-  #define BUILD_ROTATION qabc_rotation(&rotation, \
+  #define BUILD_ROTATION() qabc_rotation(&rotation, \
     theta, phi, psi, local_values.table.theta, \
     local_values.table.phi, local_values.table.psi)
-  #define APPLY_ROTATION qabc_apply(rotation, qx, qy, &qa, &qb, &qc)
-  #define CALL_KERNEL CALL_IQ_ABC(qa, qb, qc, local_values.table)
+  #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc)
+  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table)
 
 #endif
@@ -538 +561 @@
     if (weight0 > cutoff) {
       pd_norm += weight0 * CALL_VOLUME(local_values.table);
-      BUILD_ROTATION;
+      BUILD_ROTATION();
 
 #ifndef USE_OPENCL
@@ -549 +572 @@
       {
 
-        FETCH_Q;
-        APPLY_ROTATION;
+        FETCH_Q();
+        APPLY_ROTATION();
 
         // ======= COMPUTE SCATTERING ==========
         #if defined(MAGNETIC) && NUM_MAGNETIC > 0
-        // Constant across orientation, polydispersity for given qx, qy
-        double scattering = 0.0;
-        // TODO: what is the magnetic scattering at q=0
-        const double qsq = qx*qx + qy*qy;
-        if (qsq > 1.e-16) {
-          double p[4];  // dd, du, ud, uu
-          p[0] = (qy*cos_mspin + qx*sin_mspin)/qsq;
-          p[3] = -p[0];
-          p[1] = p[2] = (qy*sin_mspin - qx*cos_mspin)/qsq;
-
-          for (int index=0; index<4; index++) {
-            const double xs = spins[index];
-            if (xs > 1.e-8) {
-              const int spin_flip = (index==1) || (index==2);
-              const double pk = p[index];
-              for (int axis=0; axis<=spin_flip; axis++) {
-                #define SLD(_M_offset, _sld_offset) \
-                    local_values.vector[_sld_offset] = xs * (axis \
-                    ? (index==1 ? -values[_M_offset+2] : values[_M_offset+2]) \
-                    : mag_sld(qx, qy, pk, values[_M_offset], values[_M_offset+1], \
-                              (spin_flip ? 0.0 : values[_sld_offset+2])))
-                #define M1 NUM_PARS+5
-                #if NUM_MAGNETIC==1
-                  SLD(M1, MAGNETIC_PAR1);
-                #elif NUM_MAGNETIC==2
-                  SLD(M1, MAGNETIC_PAR1);
-                  SLD(M1+3, MAGNETIC_PAR2);
-                #elif NUM_MAGNETIC==3
-                  SLD(M1, MAGNETIC_PAR1);
-                  SLD(M1+3, MAGNETIC_PAR2);
-                  SLD(M1+6, MAGNETIC_PAR3);
-                #else
-                  for (int sk=0; sk<NUM_MAGNETIC; sk++) {
-                      SLD(M1+3*sk, slds[sk]);
-                  }
-                #endif
-                #undef SLD
-                scattering += CALL_KERNEL;
+          // Compute the scattering from the magnetic cross sections.
+          double scattering = 0.0;
+          const double qsq = qx*qx + qy*qy;
+          if (qsq > 1.e-16) {
+            // TODO: what is the magnetic scattering at q=0
+            const double px = (qy*cos_mspin + qx*sin_mspin)/qsq;
+            const double py = (qy*sin_mspin - qx*cos_mspin)/qsq;
+
+            // loop over uu, ud real, du real, dd, ud imag, du imag
+            for (int xs=0; xs<6; xs++) {
+              const double xs_weight = spins[xs];
+              if (xs_weight > 1.e-8) {
+                // Since the cross section weight is significant, set the slds
+                // to the effective slds for this cross section, call the
+                // kernel, and add according to weight.
+                for (int sk=0; sk<NUM_MAGNETIC; sk++) {
+                  const int32_t mag_index = NUM_PARS+5 + 3*sk;
+                  const int32_t sld_index = slds[sk];
+                  const double mx = values[mag_index];
+                  const double my = values[mag_index+1];
+                  const double mz = values[mag_index+2];
+                  local_values.vector[sld_index] =
+                    mag_sld(xs, qx, qy, px, py, values[sld_index+2], mx, my, mz);
+                }
+                scattering += xs_weight * CALL_KERNEL();
               }
             }
           }
-        }
         #else  // !MAGNETIC
-        const double scattering = CALL_KERNEL;
+          const double scattering = CALL_KERNEL();
         #endif // !MAGNETIC
 //printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight, weight0);
 
         #ifdef USE_OPENCL
-        this_result += weight0 * scattering;
+          this_result += weight0 * scattering;
         #else // !USE_OPENCL
-        result[q_index] += weight0 * scattering;
+          result[q_index] += weight0 * scattering;
         #endif // !USE_OPENCL
       }
@@ -627 +636 @@
   PD_CLOSE(4)
 #endif
-
-// clear the macros in preparation for the next kernel
-#undef PD_INIT
-#undef PD_OPEN
-#undef PD_CLOSE
-#undef FETCH_Q
-#undef BUILD_ROTATION
-#undef APPLY_ROTATION
-#undef CALL_KERNEL
 
 // Remember the current result and the updated norm.
@@ -646 +646 @@
 //printf("res: %g/%g\n", result[0], pd_norm);
 #endif // !USE_OPENCL
-}
+
+// clear the macros in preparation for the next kernel
+#undef PD_INIT
+#undef PD_OPEN
+#undef PD_CLOSE
+#undef FETCH_Q
+#undef BUILD_ROTATION
+#undef APPLY_ROTATION
+#undef CALL_KERNEL
+}