On Wed, 25 Feb 1998, Edwin E Evans Jr. wrote:

>  I have some expressions almost ready (on paper). I also have a copy of  
>  Greg Daues' Minimal Distortion routine from the H code.  It solves a  
>  very similar problem.   


Great! What type of formula are you using, and how are you going to 
produce the code? I would recommend using Mathematica or so.
The first thing to try might be to just lift the formula from 
Greg's code.

It would be useful to have some sort of flexibility, because for
example for maximal slicing it does matter whether one uses the "upper" or
"lower" metric form of the Laplacian.

For minimal distortion, what looks good on paper to me is to actually
finite difference Lie_beta g_ab because that's how the shift appears in
the ADM equations (with Wald's sign conventions, I think there is a sign 
difference to York somewhere), but there are many possibilities:

Laplace_l beta^a 
=
Laplace beta^a + 1/3 Del^a(Del_b beta^b) + R^a_b beta^b        (1)
=
(g^ac g^bd + g^ad g^bc - 2/3 g^ab g^cd) Del_b Del_c beta_d     (2)
=
(g^ac g^bd + g^ad g^bc - 2/3 g^ab g^cd) Del_b Lie_beta g_cd    (3)

where I have used
Lie_beta g_cd = Del_a beta_b + Del_b beta_a
You may want to check those formulas, I just did them in my notes
starting from (1).

(1) has the advantage that R_ab may be reused for the Hamiltonian.
(2) is cleanly separating taking finite differences in beta from the 
    contractions
(3) uses Lie_beta as in ADM.



>  What I need is the form you want the routine(s) to  
>  take.   

Suppose you want to solve L(u) = 0.
The basic routine needed for relaxation returns L(u) at a given point, say
indexed by ijk, AND the "diagonal derivative" for an approximation of the 
inverse Jacobian, which is trivial in the linear case. Have a look at 
stencil.c of BAM, e.g.


/* linear:
   L u  =  gi del del u  +  gg del u  +  M u  +  N
*/
void stencil_linear(struct box *g, int nfu, int ijk, 
		    double *lu, double *lii) 
{
  int di = g->di, dj = g->dj, dk = g->dk;
  double dx = g->d, dy = g->d, dz = g->d;
  double *u = g->store[nfu];

  /* inlined by hand, just to be sure = setfieldpointers_linear */
  int i = MGF*g->ncomps;
  double  *M = g->store[i++];
  double  *N = g->store[i++];
  double  *gi11 = g->store[i++];
  double  *gi12 = g->store[i++];
  double  *gi13 = g->store[i++];
  double  *gi22 = g->store[i++];
  double  *gi23 = g->store[i++];
  double  *gi33 = g->store[i++];
  double  *gg1  = g->store[i++];
  double  *gg2  = g->store[i++];
  double  *gg3  = g->store[i++];

*lu
=
M[ijk]*u[ijk] + N[ijk] + 
   gg1[ijk]*(-u[-di + ijk] + u[di + ijk])/(2*dx) + 
   gg2[ijk]*(-u[-dj + ijk] + u[dj + ijk])/(2*dy) + 
   gg3[ijk]*(-u[-dk + ijk] + u[dk + ijk])/(2*dz) + 
   gi12[ijk]*(u[-di - dj + ijk] - u[di - dj + ijk] - u[-di + dj + ijk] + 
       u[di + dj + ijk])/(2*dx*dy) + 
   gi13[ijk]*(u[-di - dk + ijk] - u[di - dk + ijk] - u[-di + dk + ijk] + 
       u[di + dk + ijk])/(2*dx*dz) + 
   gi23[ijk]*(u[-dj - dk + ijk] - u[dj - dk + ijk] - u[-dj + dk + ijk] + 
       u[dj + dk + ijk])/(2*dy*dz) + 
   gi11[ijk]*(-2*u[ijk] + u[-di + ijk] + u[di + ijk])/(dx*dx) + 
   gi22[ijk]*(-2*u[ijk] + u[-dj + ijk] + u[dj + ijk])/(dy*dy) + 
   gi33[ijk]*(-2*u[ijk] + u[-dk + ijk] + u[dk + ijk])/(dz*dz)
;

if (!lii) return;
    
*lii 
=  
M[ijk] - 2*(gi11[ijk]/(dx*dx) + gi22[ijk]/(dy*dy) + gi33[ijk]/(dz*dz))
; 

}



This should be self-explanatory: all it comes down to are finite
differences for the partial derivatives, and the coefficients have been
precomputed and stored during setup (i.e. gi is the inverse metric, gg is
a Christoffel that has been contracted once).

You should
1. try to organize the vector Laplacian in a manner such that (most of) 
   the coefficients are precomputed
2. try to avoid taking finite differences of finite differences, which 
   looses accuracy and requires two sweeps over the grid (that is, if
   we don't decide that there is a good reason to go with Lie_beta)
3. don't worry about the lii term for now, will do that later
4. compute all 3 or 4 vector components of the coupled equation during one
   call

  

>  Do you want each of the hamiltionian and momentum constraints  
>  seperately or all at once.   

Compute them together, this is the more common case and it is easier 
to put an if clause to cut out the Hamiltonian later.



>  Secondly, since the Minimal Distortion  
>  equation is so similar to the Momentum constraints I think it would be  
>  useful to add the ability to solve equations of that form generically to  
>  cactus. Much in the same way that the hamiltionian constraint and  
>  maximal slicing just interface the linear solver.  Tell me what you  
>  think.    
>  				Ed Evans 
>   


Yes. It's all about computing the coefficients. Make three routines:
for computing the Laplace_l coefficients, and for computing the sources
for constraints and gauge choices. Given these coefficients, one can call
different elliptic solvers (that support vectors).


The basic call structure should be:

compute coefficients
call standard vector elliptic solver
     which at the innermost loop calls the generic stencil operator  



-------------
All the above are just suggestions based on what worked well for 
maximal slicing, and linear and non-linear initial data equations, all
scalar. All comments are welcome.

Let me know how it is going!
Bernd

Created for the WUGRAV CoCoBoard Page Projects Page.

Created by The CoCoBoard.