Added full_chain_stiffness() to handle stiffness with piston tension.

Currently just sets a hardcoded stiffness dx, but that's pretty ugly.
A slightly better solution would be a user-controlled option, but that
is low on my agenda...

diff --git a/src/sawsim.nw b/src/sawsim.nw
index c247d7c49f48dca8b1abe9a619c55a80f2f0c546..87ed220959228790a7f1b41a20171fb1a358a64e 100644 (file)
--- a/src/sawsim.nw
+++ b/src/sawsim.nw
@@ -2316,6 +2316,7 @@ the implementation in [[tension_balance.c]], and the unit testing in
 <<x of x0>>
 <<group stiffness function>>
 <<chain stiffness function>>
+<<full chain stiffness function>>
 <<tension balance function>>
 @
 
@@ -2431,9 +2432,14 @@ double tension_balance(int num_tension_groups,
 
   F = (*tension_handlers[0])(xo, params[0]);
 
-  if (stiffness != NULL)
+  if (stiffness != NULL) {
     *stiffness = chain_stiffness(&x_xo_data, min_relx);
-
+    if (*stiffness == 0.0) { /* re-calculate based on full chain */
+      *stiffness = full_chain_stiffness(num_tension_groups, tension_handlers,
+                                        inverse_tension_handlers, params,
+                                        last_x, x, min_relx, F);
+    }
+  }
   return F;
 }
 
@@ -2499,6 +2505,10 @@ double x_of_xo(double xo, void *pdata)
 }
 @
 
+Determine the stiffness of a single tension group by taking a small
+step [[dx]] back from the position [[x]] for which we want the
+stiffness.  The touchy part is determining a good [[dx]] and ensuring
+the whole interval is on [[x>=0]].
 <<group stiffness function>>=
 double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx)
 {
@@ -2524,10 +2534,19 @@ double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double
 }
 @
 
+Attempt to calculate the chain stiffness by adding group stiffnesses
+as springs in series.  In the case where a group has undefined
+stiffness (e.g. piston tension, Section \ref{sec.piston_tension}),
+this algorithm breaks down.  In that case, [[tension_balance()]]
+(\ref{sec.tension_balance}) should use [[full_chain_stiffness()]]
+which uses the full chain's $F(x)$ rather than that of the individual
+domains'.  We attempt the series approach first because, lacking the
+numerical inversion inside [[tension_balance()]], it is both faster
+and more accurate than the full-chain derivative.
 <<chain stiffness function>>=
 double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
 {
-  double stiffness;
+  double stiffness, gstiff;
   int i;
   /* add group stiffnesses in series */
   for (i=0; i < x_xo_data->n_groups; i++) {
@@ -2535,7 +2554,10 @@ double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
     fprintf(stderr, "%s:\tgetting stiffness of active state %d of %d for x[%d]=%g, relx=%g\n", __FUNCTION__, i, x_xo_data->n_groups, i, x_xo_data->xi[i], relx);
     fflush(stderr);
 #endif
-    stiffness += 1.0/group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
+    gstiff = group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
+    if (gstiff == 0.0);
+      return 0.0;
+    stiffness += 1.0/gstiff;
   }
   stiffness = 1.0 / stiffness;
   return stiffness;
@@ -2543,6 +2565,63 @@ double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
 
 @
 
+Determine the chain stiffness using only [[tension_balance()]].  This
+works around difficulties with tension models that have undefined
+stiffness (see discussion for [[chain_stiffness()]] above).
+<<full chain stiffness function>>=
+double full_chain_stiffness(int num_tension_groups,
+                       one_dim_fn_t **tension_handlers,
+                      one_dim_fn_t **inverse_tension_handlers,
+                       void **params, /* array of pointers to tension_handler_data_t */
+                      double last_x,
+                      double x,
+                      double relx,
+                      double F /* already evaluated F(x) */)
+{
+  double dx, xl, Fl, stiffness;
+
+  assert(x >= 0);
+  assert(relx > 0 && relx < 1);
+
+  /* Other option for dx that we ignore because of our poor tension_balance()
+   * resolution (i.e. bad slopes if you zoom in too close):
+   *     if (last_x != -1.0 && last_x != x)  // last_x limited
+   *       dx fabs(x-last_x);
+   *     else
+   *       dx = HUGE_VAL;
+   *     if (1==1) {                 // constant max-value limited
+   *       dx_p = 1e-12;
+   *       dx = MIN(dx, dx_p);
+   *     }
+   *     if (x != 0 && relx != 0) {  // relx limited
+   *       dx_p = x*relx;
+   *       dx = MIN(dx, dx_p);
+   *     }
+   * TODO, determine which of (min_relx, min_rely, max_steps) is actually
+   * limiting tension_balance.
+   */
+  dx = 1e-12; /* HACK, how to get length scale? */
+
+  xl = x-dx;
+  if (xl >= 0) {
+    Fl = tension_balance(num_tension_groups, tension_handlers,
+                         inverse_tension_handlers, params,
+                         last_x, xl, NULL);
+  } else {
+    xl = x;
+    Fl = F;
+    x += dx;
+    F = tension_balance(num_tension_groups, tension_handlers,
+                       inverse_tension_handlers, params,
+                       last_x, x, NULL);
+  }
+  stiffness = (F-Fl)/dx;
+  assert(stiffness >= 0);
+  return stiffness;
+}
+
+@
+
 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
 number of steps for monotonically increasing $f(x)$.  Simple
 bisection, so it's robust and converges fairly quickly.  You can set
@@ -4361,6 +4440,12 @@ tension at $x=L$ is undefined, but may be any positive value.  The
 model is called the ``piston'' model because it resembles a
 frictionless piston sliding in a rigid cylinder of length $L$.
 
+Note that because of it's infinte stiffness at $x=L$, fully extended
+domains with this tension model will be identical to completely rigid
+domains.  However, there is a distinction between this model and rigid
+domains for $x<L$, because any reactions that occur at $F=0$
+(e.g. refolding) will have more time to occur.
+
 Because the tension at $x=L$ is undefined, the user must make sure
 domains with this tension model are never the initial active group,
 because it would break [[tension_balance()]] in [[find_tension()]]