[CIG-MC] CitcomCU stress dependent rheology

[Eh Tan]

On 3/19/10 8:29 AM, Ethan Fahy wrote:
> Hello Eh,
>
> In the CitcomCU readme it says to be careful using stress dependent
> rheology, and indeed when I turn it on and set the stress exponents to
> anything other than 1.0 lots of nans come up in the screen output.  I
> have seen in multiple papers the use of a lower activation energy to
> mimic the effects of stress dependent rheology.  Is true stress
> dependent rheology in CitcomCU something that has been done
> successfully?  If so are there any hints you might be able to give me
> to get it working correctly?  Thanks very much.
>
> -Ethan
>    

Hi Ethan,
(cc to cig-mc)

The non-newtonian solver in CitcomCU-1.0.3 has a bug. Magali Billen 
reported the bug to me off-list in November. The bug was fixed in the 
repository, but I forgot to announce it to the mailing list. I am sorry 
about it.

Here is the fix of the bug in diff format.

--- Drive_solvers.c     (revision 16060)
+++ Drive_solvers.c     (revision 16061)
@@ -77,15 +77,15 @@
               damp = 0;
             }
           } /* end powerlaw */
-               oldU = (double *)malloc((neq + 2) * sizeof(double));
-               for(i = 1; i <= neq; i++)
+               oldU = (double *)malloc(neq * sizeof(double));
+               for(i = 0; i < neq; i++)
                         oldU[i] = 0.0;
                 visits++;
         }

         dUdot_mag = 0.0;

-       delta_U = (double *)malloc((neq + 2) * sizeof(double));
+       delta_U = (double *)malloc(neq * sizeof(double));

         /* FIRST store the old velocity field */

@@ -124,10 +124,10 @@

                   if(damp){
                     /* add some of the old solution */
-                   for(i = 1; i <= neq; i++)
+                   for(i = 0; i < neq; i++)
                       E->U[i] = alpha * E->U[i] + alpha1 * oldU[i];
                   }
-                       for(i = 1; i <= neq; i++)
+                       for(i = 0; i < neq; i++)
                         {
                                 delta_U[i] = E->U[i] - oldU[i];
                                 oldU[i] = E->U[i];

-- 
Eh Tan
-------
Postdoctoral Fellow
Institute of Geophysics
University of Texas, Austin
J.J. Pickle Research Campus, Bldg 196 (ROC)
10100 Burnet Rd. (R2200),
Austin, TX 78758-4445
(512) 475-9568