SchwarzschildF5.c

Writing a time-dependent metric tensor field on a uniform grid.

 
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <hdf5.h>
 
#include <F5/F5X.h>
#include <F5/F5uniform.h>
 
/* #include <malloc.h> */
 
int     main(int argc, char*argv[])
{
        double  time = 0.0;
        hid_t   file_id; 
 
        hsize_t dims[3] = {10,20,30};
        /*
        Note: do not use multidimensional C arrays.
        In the C notation, an array data[X][Y][Z] implies that
        the Z coordinate runs fastest and succeeding points in z are
        close in memory, whereas points close in X are separated by slices.
        Using C arrays, it is not possible to perform XY operations trivially.
        Thus in F5 we prefer to use FORTRAN order, which in C looks like
        data[Z][Y][X] , or, better, use one-dimensional arrays and do
        indexing yourself. This the more real-world case anyway, because
        the dimension of a multidimensional array must be known at compile-time,
        which is normally not the case.
        */
        F5_metric3_float_t m[10*20*30]; 
        F5_vec3_point_t origin; 
        F5_vec3_float_t delta; 
 
        F5Path*myPath;
                
        unlink("metric.f5");
        file_id = F5append( "metric.f5" ); 
 
        origin.x = -1;               origin.y = -1;                     origin.z = -1; 
        delta .x = 2.0/(int)dims[0]; delta .y = 2.0/(int)dims[1];       delta .z = 2.0/(int)dims[2];    
        
        for(time=0.0; time<10; time+=0.5)
        {
                int     i,j,k;          
                double  M = (10 - time)/10; 
                F5_metric3_float_t *g = m;
 
                for(k=0; k<dims[2]; k++)
                for(j=0; j<dims[1]; j++)
                for(i=0; i<dims[0]; i++)
                {
                        float x, y, z, r, psi;
                
                        x = origin.x + i * delta.x; 
                        y = origin.y + j * delta.y; 
                        z = origin.z + k * delta.z; 
                        r = sqrt(x*x+y*y+z*z); 
                        
                        psi = 1 + M / (2*r); 
                        psi *= psi;     // psi^2 
                        psi *= psi;     // psi^4 
 
        
 
                        g->gxx = psi;
                        g->gyy = psi; 
                        g->gzz = psi; 
                        g->gxy = 0; 
                        g->gyz = 0; 
                        g->gxz = 0; 
                        
                        g++;
                } 
 
                myPath = F5Fwrite_uniform_cartesian3D(file_id, time,
                                                      "schwarzschild",
                                                      &origin, &delta,
                                                      dims,
                                                      "metric",
                                                      F5T_METRIC33_FLOAT,
                                                      m, 0, F5P_DEFAULT); 
                F5close(myPath);
        }
        H5Fclose(file_id);      
        return 0;
}