[CIG-MC] Ridge rots and blow ups in CitcomS
Shijie Zhong
Shijie.Zhong at Colorado.Edu
Sat Oct 13 09:30:48 PDT 2007
Hi Thorsten,
Glad to see that the iteration problem is resolved so quickly.
For benchmark results for CitcomS, you can already check the following website for ~20 different cases with highly variable viscosity:
http://www.geodynamics.org/cig/workinggroups/mc/workarea/benchmark/3dconvention/
Cheers,
Shijie
Shijie Zhong
Department of Physics
University of Colorado at Boulder
Boulder, CO 80309
Tel: 303-735-5095; Fax: 303-492-7935
Web: http://anquetil.colorado.edu/szhong
---- Original message ----
>Date: Sat, 13 Oct 2007 08:37:43 -0700
>From: Thorsten Becker <twb at usc.edu>
>Subject: Re: [CIG-MC] Ridge rots and blow ups in CitcomS
>To: Shijie Zhong <Shijie.Zhong at Colorado.Edu>
>Cc: Eh Tan <tan2 at geodynamics.org>,cig-mc at geodynamics.org
>
>Dear Shijie,
>
>Thanks for the explanation, great to hear that there's a benchmark
>going on. I am working on the SVN version of CitcomS when adding my
>modifications, that's why I got the update already.
>
>As you know, following your 2001 EPSL paper, I'm a big fan of the net
>rotations myself, and noticed that those can be quite strong
>particularly for plastic yielding spherical models, quite interesting.
>That's why I put in the lib/Determine_net_rotation.c a while back.
>
>However, as to be expected, your implementation of the net rotation
>removal should be much quicker than mine. Also, I forgot about the time
>step issue and agree that it's much better to remove the net rotation
>after every velocity solution, rather than removing only at output.
>
>The problem that caused the blowup when your routine was called within
>CIG CitcomS was very simple; I would fix it in the SVN version but
>cannot access my cluster right now:
>
>VV should be declared float rather than double in remove_rigid_rot
>
>float VV[4][9]; (l. 819 of Global_operations.c)
>
>because (at least in my version)
>
>
>void velo_from_element(E,VV,m,el,sphere_key)
> struct All_variables *E;
> float VV[4][9];
> int el,m,sphere_key;
>
>in Nodal_mesh.c
>
>It would be very easy to allow the compiler to catch those declarative
>mismatches from now on for eternity, if we would switch the function
>declarations from old-school type to the
>
>void velo_from_element(struct All_variables *,
> float [4][9],int ,int ,int);
>
>type. Each file could include a "functions.h" file which has all the
>function declarations. This could be generated automatically by
>
> cproto $(INCLUDES) -f2 -q *.c | grep -v "void main(" | grep -v "int
>main(" > functions.h
>
>In fact, Luiz once made such changes, but then they must have been taken
>out of the major source tree of CitcomS for some reason.
>
>Cheers
>
>T
>
>
>
>
>
>
>On Sat, 2007-10-13 at 03:18 -0600, Shijie Zhong wrote:
>> Thorsten,
>>
>> remove_rigid_rot() was added only ~ a week ago (you really keep things updated). In the new release (possibly next week?), the purpose of this routine will be explained. I am working on a benchmark paper (would have been done already if not for this rigid body rotation) that would detail this issue. I can explain to you (and those who may be interested) briefly here.
>>
>> For the stokes' flow in a spherical shell, it is possible to add an arbitrary ridge body rotation to the flow and the governing equations will still be satisfied (i.e., the rigid body rotation is unconstrained). In most calculations (as far as I know), this mode of rigid body rotation is never important. However, for some reasons, it may show up in some other cases (e.g., thermochemical convection, as Allen McNamara discovered, or mobile-lid convection with continental keels). When the rigid body rotation is present in a calculation, it can seriously and unnecessarily reduce time step (delta t) in convection calculation, as delta t is determined by maximum velocity, and it also leads to ambiguity in velocities (i.e., reference frame issue, important for plate tectonics related topics and benchmark purposes).
>>
>> I think that a proper reference frame is no-net-rotation for the whole MANTLE (different from no-net-rotation for lithosphere as used sometimes in plate tectonics related studies). In the past, as what you have done, I have used a post-processing routine to remove the rigid body rotation for the whole mantle (e.g., Zhong [2001]). Allen then found that the rigid body rotation sometimes gave him too small delta t in time-dependent thermochemical convection calculations, so he coded a routine into CitcomS (thermochemical one in 2004) that can be called every time step.
>>
>> I recently "re-learnt" the lesson when working on a benchmark paper (more complete description on this issue is in that paper). I "wasted" more than a month in trying to figure out why two resolution test cases for thermochemical convection that differ only in resolution have very different Vrms, while heat flux is nearly the same. It turns out that the high resolution case somehow has some rigid body rotation.
>> At the end, I think that we should include this routine and that the rigid body rotation should be removed every time step, although it does not affect the dynamics. This also gives us a reference frame for all the velocities (i.e., no-net-rotation for the mantle). There should be two routines in CitcomS right now, one is Allen's original one, and the other (being called) is a simplified version from Allen's.
>>
>> Now to the non-Newtonian issue, I cannot see why this routine should not mess up the non-Newtonian iterations. There might be some incompartibility issue between this routine and CIG's version (?). Perhaps, you can turn it off for now.
>>
>>
>> Shijie
>>
>>
>> Shijie Zhong
>> Department of Physics
>> University of Colorado at Boulder
>> Boulder, CO 80309
>> Tel: 303-735-5095; Fax: 303-492-7935
>> Web: http://anquetil.colorado.edu/szhong
>>
>>
>> ---- Original message ----
>> >Date: Fri, 12 Oct 2007 17:59:15 -0700
>> >From: Thorsten Becker <twb at usc.edu>
>> >Subject: [CIG-MC] Ridge rots and blow ups in CitcomS
>> >To: Eh Tan <tan2 at geodynamics.org>
>> >Cc: cig-mc at geodynamics.org
>> >
>> >Hi Eh,
>> >
>> >I am a bit confused about the net rotation removal routine
>> >
>> >remove_rigid_rot in Global_operations.c
>> >
>> >in CitcomS. Has this always been there, or is this a relatively new
>> >addition?
>> >
>> >I had recently also built in a (slightly different) net rotation
>> >removal routine as
>> >
>> >lib/Determine_net_rotation.c:53:double
>> >determine_model_net_rotation(struct All_variables *,double *);
>> >
>> >which I activate only in the gzipped output routine in a post-processing
>> >step if
>> >
>> >gzdir_rnr=on
>> >
>> >
>> >and not at every time step. I noticed the new routine because it lead
>> >to a blow up of velocities if stress-dependent viscosity iterations
>> >were performed after a tic_method=-1 restart, for reasons unclear to me
>> >right now.
>> >
>> >Could we maybe set
>> >
>> >remove_rigid_rotation
>> >
>> >to an "off" default for now?
>> >
>> >Thanks for any comments
>> >
>> >T
>> >
>> >
>> >
>> >
>> >--
>> > Thorsten W Becker Department of Earth Sciences
>> > University of Southern California p: 213.740.8365
>> > Los Angeles CA 90089-0740 http://geodynamics.usc.edu/~becker
>> >
>> >_______________________________________________
>> >CIG-MC mailing list
>> >CIG-MC at geodynamics.org
>> >http://geodynamics.org/cgi-bin/mailman/listinfo/cig-mc
>--
> Thorsten W Becker Department of Earth Sciences
> University of Southern California p: 213.740.8365
> Los Angeles CA 90089-0740 http://geodynamics.usc.edu/~becker
>
More information about the CIG-MC
mailing list