[CIG-MC] CitcomS 3.1.1 Cookbook6 not converging

[Robert Moucha]

Thanks Eh,

I will give this a try.

Cheers
Rob


On Mon, Sep 14, 2009 at 7:49 PM, Eh Tan <tan2 at geodynamics.org> wrote:
> Hi Robert,
>
> There is a bug in one pseudo free surface subroutine in
> lib/Element_calculations.c. Please replace the subroutine
> assemble_forces_pseudo_surf() with the code attached below. This should
> fix the problem.
>
>
> Cheers,
> Eh
>
> ------------------------------------
>
> void assemble_forces_pseudo_surf(E,penalty)
>     struct All_variables *E;
>     int penalty;
> {
>  double elt_f[24];
>  int m,a,e,i;
>
>  void get_buoyancy();
>  void get_elt_f();
>  void get_elt_tr_pseudo_surf();
>  void strip_bcs_from_residual();
>  double global_vdot();
>
>  const int neq=E->lmesh.neq;
>  const int nel=E->lmesh.nel;
>  const int lev=E->mesh.levmax;
>
>  get_buoyancy(E,E->buoyancy);
>
>  for(m=1;m<=E->sphere.caps_per_proc;m++)    {
>
>    for(a=0;a<neq;a++)
>      E->F[m][a] = 0.0;
>
>    for (e=1;e<=nel;e++)  {
>      get_elt_f(E,e,elt_f,1,m);
>      add_force(E, e, elt_f, m);
>    }
>
>    /* for traction bc */
>    for(i=1; i<=E->boundary.nel; i++) {
>      e = E->boundary.element[m][i];
>
>      for(a=0;a<24;a++) elt_f[a] = 0.0;
>      for(a=SIDE_BEGIN; a<=SIDE_END; a++)
>    get_elt_tr_pseudo_surf(E, i, a, elt_f, m);
>
>      add_force(E, e, elt_f, m);
>    }
>  }       /* end for m */
>
>  (E->solver.exchange_id_d)(E, E->F, lev);
>  strip_bcs_from_residual(E,E->F,lev);
>
>  /* compute the norm of E->F */
>  E->monitor.fdotf = sqrt(global_vdot(E, E->F, E->F, lev));
>
>  if(E->parallel.me==0) {
>      fprintf(stderr, "Momentum equation force %.9e\n",
>              E->monitor.fdotf);
>      fprintf(E->fp, "Momentum equation force %.9e\n",
>              E->monitor.fdotf);
>  }
>
>  return;
> }
>
>
>
> Robert Moucha wrote:
>> Hi All,
>>
>> I was wondering if anyone has ran the Cookbook6 (pseudo surface)
>> example using 3.1.1 with the default cookbook6.cfg settings. All runs
>> well for the free slip case, but for the pseudo surface case I get the
>> following output after an hour and a half.
>>
>> Any ideas? I have not played with the settings because I'm under the
>> impression that it should run straight out of "the box".
>> Thanks
>>
>> Rob Moucha
>>
>>
>> Problem has 61 x 25 x 61 nodes per cap, 93025 nodes and 86400 elements in total
>>  *  residual (100000) = 8.112e-161 from 8.564e+02 to 0.000e+00 in 4719.55 secs
>> Warning: solver not converging! 0
>> (000) 4719.6 s v=3.530406e+03 p=0.000000e+00 div/v=5.19e+00
>> dv/v=1.00e+00 dp/p=1.00e+00 step 0
>>  *  residual (100000) = 7.674e-161 from 1.393e+04 to 0.000e+00 in 4621.91 secs
>> Warning: solver not converging! 1
>> (001) 9341.5 s v=2.915688e+03 p=2.625429e+05 div/v=2.20e+00
>> dv/v=3.28e-01 dp/p=1.00e+00 step 0
>>  *  residual (100000) = 7.298e-161 from 7.563e+03 to 0.000e+00 in 4686.65 secs
>> Warning: solver not converging! 1
>> (002) 14028.2 s v=2.918320e+03 p=2.538605e+05 div/v=8.71e-01
>> dv/v=8.37e-02 dp/p=4.00e-01 step 0
>> _______________________________________________
>> CIG-MC mailing list
>> CIG-MC at geodynamics.org
>> http://geodynamics.org/cgi-bin/mailman/listinfo/cig-mc
>>
>
> --
> Eh Tan
> Staff Scientist
> Computational Infrastructure for Geodynamics
> California Institute of Technology, 158-79
> Pasadena, CA 91125
> (626) 395-1693
> http://www.geodynamics.org
>
>



-- 
GEOTOP - Département des Sciences de la Terre et de l'Atmosphère
Université du Québec à Montréal
CP 8888, succursale Centre-Ville
Montréall, Québec
Canada  H3C 3P8
Tel:     (1-514) 987-3000, ext 1554#
FAX:     (1-514) 987-3635