Turbulence.c File Reference

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

Include dependency graph for Turbulence.c:

Go to the source code of this file.

Macros
#define	Points_X   128
#define	Points_Y   128
#define	x_end   1.
#define	Norm   .1 /* Normalization of the spectrum power-law */
#define	Spectral_Index   -5./3. /* Index of the spectrum power-law: Kolmogorov -> 5/3 */
#define	IR_Cutoff   2. /* Needed to avoid an IR divergence; 2 is a good value */
#define	Minimum_Lambda   4 /* Grid points corresponding to minimum wavelength */
#define	frand()   ((double)(rand()+1)/(RAND_MAX+1.0))

Functions
double **	array_2D (int nx, int ny)
int	main (void)

Macro Definition Documentation

#define frand

(

)

((double)(rand()+1)/(RAND_MAX+1.0))

Definition at line 20 of file Turbulence.c.

#define IR_Cutoff   2. /* Needed to avoid an IR divergence; 2 is a good value */

Definition at line 17 of file Turbulence.c.

#define Minimum_Lambda   4 /* Grid points corresponding to minimum wavelength */

Definition at line 18 of file Turbulence.c.

#define Norm   .1 /* Normalization of the spectrum power-law */

Definition at line 15 of file Turbulence.c.

#define Points_X   128

Definition at line 12 of file Turbulence.c.

#define Points_Y   128

Definition at line 13 of file Turbulence.c.

#define Spectral_Index   -5./3. /* Index of the spectrum power-law: Kolmogorov -> 5/3 */

Definition at line 16 of file Turbulence.c.

#define x_end   1.

Definition at line 14 of file Turbulence.c.

Function Documentation

double ** array_2D	(	int	nx,
		int	ny
	)

Definition at line 111 of file Turbulence.c.

    {
    int i;
    double **m;
    
    m = (double **)malloc((size_t)nx*sizeof(double *));
    m[0] = (double *)malloc((size_t)nx*ny*sizeof(double));

    for(i = 1; i < nx; i++) m[i] = m[i - 1] + ny;
    
    return m;
    }

Here is the caller graph for this function:

int main

(

void

)

Definition at line 24 of file Turbulence.c.

{
  int i, j, k, i2, j2, k2, te0, NXLim, NYLim;
  double y_end, xdbl, NXLim_2, NYLim_2;
  double twoPI, x1r, x2r, twopi_a, twopi_b, twopi_c, scrh1, scrh2;
  double xL, xR;
  double **Amplitude, **Phase;
    
  FILE *Out;

/* --------------------------------------------------------
             Grid generation
   -------------------------------------------------------- */

  Out = fopen("grid0.out", "w");
  fprintf(Out, "# GEOMETRY:   CARTESIAN\n");

  fprintf(Out, "%d\n",Points_X);
  for(i2 = 0; i2 < Points_X; i2++){
    xL = -0.5 + (double)(x_end*(i2-0.5)/(Points_X - 1.));
    xR = -0.5 + (double)(x_end*(i2+0.5)/(Points_X - 1.));
    fprintf (Out, "%d   %12.6e  %12.6e \n",i2+1, xL, xR);
  }

  fprintf(Out, "%d\n",Points_Y);
  for(j2 = 0; j2 < Points_Y; j2++){
    xL = -0.5 + (double)(x_end*(j2-0.5)/(Points_Y - 1.));
    xR = -0.5 + (double)(x_end*(j2+0.5)/(Points_Y - 1.));
    fprintf (Out, "%d   %12.6e  %12.6e \n",j2+1, xL, xR);
  }

  fprintf(Out, "%d\n",1);  /* write 3rd coordinates as well */
  fprintf (Out, "1  0.0  1.0\n");
  fclose(Out);

/* --------------------------------------------------------
               Data generation
   -------------------------------------------------------- */

  te0 = time(NULL);
  srand(te0);
  NXLim = 2*((int)(Points_X/2./Minimum_Lambda));
  NYLim = 2*((int)(Points_Y/2./Minimum_Lambda));    
  NXLim_2   = NXLim/2.;
  NYLim_2   = NYLim/2.;
  y_end = x_end*(Points_Y - 1.)/(Points_X - 1.);

  Amplitude = array_2D(NXLim + 1, NYLim + 1);
  Phase     = array_2D(NXLim + 1, NYLim + 1);
  twoPI     = 2.*3.14159265359;
    
  for(i = 0; i <= NXLim; i++) for(j = 0; j <= NYLim; j++){
    x1r = frand();
    x2r = frand();
    xdbl = sqrt(-2.*log(x1r))*cos(twoPI*x2r);
    Phase[i][j] = twoPI*frand();
    scrh1       = IR_Cutoff*IR_Cutoff + pow(twoPI*(i - NXLim_2)/x_end,2.)
                                      + pow(twoPI*(j - NYLim_2)/y_end,2.);
    Amplitude[i][j] = sqrt(Norm*pow(scrh1, Spectral_Index/2.))*xdbl;
  }
  Out = fopen("rho0.dbl", "wb");

  for(j2 = 0; j2 < Points_Y; j2++){
    twopi_b = twoPI/(Points_Y - 1.)*j2;
    for(i2 = 0; i2 < Points_X; i2++){
      twopi_a = twoPI/(Points_X - 1.)*i2;
      xdbl = 0.;
      for(i = NXLim + 1; i--; ){
        scrh1 = twopi_a*(i - NXLim_2); 
        for(j = NYLim + 1; j--; ){
          scrh2 = twopi_b*(j - NYLim_2);
          xdbl += Amplitude[i][j]*cos(scrh1 + scrh2 + Phase[i][j]);
        }
      }
      xdbl = exp(xdbl);
      fwrite (&xdbl, sizeof(double), 1, Out);
    }
  }

  fclose(Out);
  return;
}

Here is the call graph for this function: