[CIG-MC] CitcomS 3.1.1 Cookbook6 not converging

[Eh Tan]

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