boussinesqc.c

Test Boussinesq Linearized
Problem
Equations
- Main
- Run
- Results
- Links
- biblio

Test Boussinesq Linearized

This example is not Basilisk but standard C for M1 course. We use explicit resolution with centered finite derivatives to solve the Boussinesq linearized wave problem.

Problem

The shallow water equations are non linear bit they are not dispersives. But water waves in moderate depths are dispersives. To fix this problem, a simple way is to consider Boussinesq correction.

Equations

We look at the solution of linearized Boussinesq equations as then come from analysis (see M2 course) \displaystyle \frac{\partial u}{\partial t} = -\frac{\partial \eta }{\partial x} \displaystyle \frac{\partial \eta }{\partial t} = - \frac{\partial u}{\partial x} + \frac{1}{3} \frac{\partial^3 \eta}{\partial x^2 \partial t } note that a better point of view, which presserves mass conservation, is to write instead \displaystyle \frac{\partial \eta}{\partial t} = -\frac{\partial u }{\partial x} \displaystyle \frac{\partial u }{\partial t} = - \frac{\partial \eta}{\partial x} + \frac{1}{3} \frac{\partial^3 u}{\partial x^2 \partial t } as the problem is linearized, it is the same problem. The equations are solved using a factorisation of the \frac{\partial^2 }{\partial x^2 } term: \displaystyle \frac{\partial u}{\partial t} = -\frac{\partial \eta }{\partial t} \displaystyle (1 - \frac{1}{3} \frac{\partial^2 }{\partial x^2 } ) \frac{\partial \eta }{\partial t} = - \frac{\partial u}{\partial x} In version of (1 - \frac{1}{3} \frac{\partial^2 }{\partial x^2 } ) will be implicit whereas the scheme is explicit. Again, this is done in simple C for comparison.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define pi 3.151592

starting form here, it can be in a invSOR.h

routine alloc and free de OOura pour créeer matrice n x n

double **alloc_2d_double(int n1, int n2)
{
    double **dd, *d;
    int j; 
    dd = (double **) malloc(sizeof(double *) * n1); 
    d = (double *) malloc(sizeof(double) * n1 * n2); 
    dd[0] = d;
    for (j = 1; j < n1; j++) {
        dd[j] = dd[j - 1] + n2;
    }
    return dd;
}
void free_2d_double(double **dd)
{
    free(dd[0]);
    free(dd);
}

solve A.X = B by “JOR/SOR Jacobi and Gauss-Seidel over-relaxation Method” Alfio Quarteroni Riccardo Sacco Fausto Saleri, page 136

void invSOR(int sor,double **A,double *X,double *B,int n)
{ 
  int iter=0,i,j;
  double s,omega=.5;
  double  r0,err;
  static double *Xbiel=NULL,*R=NULL;  
  Xbiel=(double*)calloc(n,sizeof(double));
  R=(double*)calloc(n,sizeof(double)); 
  s=0;
  for(i=0;i<n;i++){
      for(j=0;j<n;j++)s+=A[i][j]*X[j];
      R[i]=B[i]-s;
   } 
   r0=0;
   for(i=0;i<n;i++)r0+=(R[i]*R[i]);
   r0=sqrt(r0);
   err=r0;
   do{  
    iter++;
    for(i=0;i<n;i++)Xbiel[i]=X[i];
    for(i=0;i<n;i++){
      s=0;  
      if(sor==1){  
      if(i>0){
      for(j=0  ;j<=i-1;j++){s+=A[i][j]*X[j]; }
      }
      if(i<n-1){
        for(j=i+1;j<n  ;j++){s+=A[i][j]*Xbiel[j];}
      }
      } 
      else{     
      for(j=0 ;j<n;j++){if(i!=j) {s+=A[i][j]*X[j]; }}
      }
      X[i]=  omega*(B[i]-s)/A[i][i]+(1-omega)*Xbiel[i] ;
     }
    for(i=0;i<n;i++){
          s=0;
      for(j=0 ;j<n;j++){  s+=A[i][j]*X[j];  }
      R[i]=B[i]-s;      
    } 
   err=0;
   for(i=0;i<n;i++)err+=(R[i]*R[i]);
    err=sqrt(err)/r0;   
      }
    while((err > 1e-12)&&(iter<20000));
  if((err > 1e-12)||(iter>20000-1)) fprintf(stderr," iter=%d err=%lf\n", iter,err*1e+12);
  free(Xbiel);
  free(R);
}

fin de la routine d’inversion SOR

end of invSOR.h

Main

Main program, value of paerameters, initial wave

int main (int argc, const char *argv[]) {
    int  i,nx,it=0;//,sor=1;
  //  char file;    
    double*x=NULL,*u=NULL,*F=NULL,*uo=NULL,*eta=NULL,*etao=NULL,*detadt=NULL;
    double dt,dx,x0,L,s,t,sigma,tmax; 
    double **a=NULL;
    FILE *g;
 
     nx = 375;
     L = 125;
     x0=-L/3.; 
     dx = L/nx;
     dt = fmin(dx,.25) ;// OK  0.1; 
     t = 0;
     tmax=60;
     sigma=1./3; //   1./3 for real Boussinesq
      
      x=     (double*)calloc(nx+1,sizeof(double));
      u=     (double*)calloc(nx+1,sizeof(double));
      eta=   (double*)calloc(nx+1,sizeof(double));
      uo=    (double*)calloc(nx+1,sizeof(double));
      etao=  (double*)calloc(nx+1,sizeof(double));
      detadt=(double*)calloc(nx+1,sizeof(double));
      F=     (double*)calloc(nx+1,sizeof(double)); 
     fprintf(stderr,"solution de  du/dt = - dh/dx; \n");
     fprintf(stderr,"solution de  dh/dt = - du/dx + sigma d^3 h/dx^2dt \n");
     fprintf(stderr,"dx=%lf,  sigma*dt//dx^3=%lf  dt/dx=%lf\n",dx,sigma*dt/dx*dx*dx,dt/dx);
// init 
     for(i=0;i<=nx;i++)
     { u[i]=0;
       x[i]= x0 +i*dx;
       uo[i]=0;
       u[i]=uo[i];
       etao[i]=exp(-pi*(x[i])*(x[i]));;
      // etao[i]=.01*(1-tanh(x[i]));;         
       u[i]=etao[i];
       uo[i]=u[i];
       eta[i]=etao[i];
     } 
     g = fopen("sol.OUT", "w");
     fclose(g);

Begin time loop while time less than tmax

  do{  
    t=t+dt; 
    it++;

Explicit resolution with centered finite derivatives. Compute \dfrac{\partial \eta}{\partial x} \sim \dfrac{ \eta(x+\Delta x)-\eta(x-\Delta x)}{2 \Delta x } so that \displaystyle u (x, t +\Delta t ) = u (x, t ) - \Delta t \dfrac{ \eta(x+\Delta x)-\eta(x-\Delta x)}{2 \Delta x }

    for(i=1;i<nx-1;i++) u[i] = uo[i]   - dt*(eta[i+1]-eta[i-1])/2/dx ;                            
            u[0]=uo[1];
            u[nx-1]=0;
            u[nx-2]=0;
// swap               
        for(i=0;i<nx;i++) uo[i]=u[i];        
// prepare    (1 - (1/3) d^2/dx^2) dh/dt = - du/dx       
      a = alloc_2d_double(nx,nx);  
      s=sigma/(dx*dx);
      for(i=2;i<nx-2;i++){a[i][i-2]=0;a[i][i-1]=-s;a[i][i]=1+2*s;a[i][i+1]=-s;a[i][i+2]=0;}
      a[0][0]=1;  
      a[1][1]=1; 
      a[nx-2][nx-2]=1;
      a[nx-1][nx-1]=1;
      detadt[0]=0;
      detadt[1]=0;
      detadt[nx-2]=0.;
      detadt[nx-1]=0.;

compute implicitely the increase in time of \eta \displaystyle F= -\dfrac{\partial u}{\partial x} \sim -\dfrac{ u(x+\Delta x)-u(x-\Delta x)}{2 \Delta x }

    for(i=1;i<nx-1;i++) F[i] = - (uo[i+1] -uo[i-1] )/dx/2;
    F[0]=0;
    F[nx-1]=0;

by inversion of (1 - \frac{1}{3} \frac{\partial^2 }{\partial x^2 } ) f = F

#ifdef bla
      invSOR(sor,a,detadt,F,nx);
    //  for(i=1;i<nx;i++)detadt[i]=F[i];
#else
      double eps=.0,detadtn=0.,omega=.25,s=1./3;
      do
      { eps=0;
           for(i=1;i<nx;i++)
          {
              detadtn = (dx*dx*F[i] + s*(detadt[i-1] + detadt[i+1]))/(2*s + dx*dx);
              eps = eps +  (detadt[i]-detadtn)*(detadt[i]-detadtn);
              detadt[i] =   omega * detadtn + (1-omega) * detadt[i] ;
          }
          // fprintf (stderr, "t=%g eps=%g \n",t, sqrt(eps));
      }while(sqrt(eps)>1.e-11);
   
#endif
      free_2d_double(a);

we update \eta (x, t+\Delta t) = \eta (x, t) + (\Delta t) f

    for(i=1;i<nx;i++) eta[i] = etao[i]   +  detadt[i]*dt ;

swap an simple BC

      eta[0] = eta[1];  
      eta[nx-1]=0;  
      eta[nx-2]=0;
      eta[nx-3]=0;
      for(i=1;i<nx;i++) etao[i]=eta[i];

affichage

    if((it%4)==0){ 
    g = fopen("sol.OUT", "a");
  for (i=0; i<=nx;i++)
  {      fprintf(g,"%lf %lf %lf \n",x[i],u[i],t);}
  fprintf(g,"\n\n");
  fclose(g);
  
  fprintf(stderr,"  t=%lf\n",t);
  
  printf(" p[%lf:%lf][%lf:%lf]  '-' u 1:2 t'u'  w l linec 3 , '-' u 1:2 t'h' w l\n ",
         x[1],x[nx],-1.,1.); 
    for (i=0; i<=nx;i++)
    {
    printf("%lf %lf \n",x[i],u[i]);}
    printf("e \n");
     for (i=0; i<=nx;i++)
    {
    printf("%lf %lf \n",x[i],eta[i]);}
    printf("e \n");
    }
  
  }while(t<tmax);

  // liberation finale
   free(x);
   free(u);
   free(uo);
   free(eta);
   free(etao);
   free(F);
   free(detadt);  
}

Run

to compile with standard C

 cc boussinesqc.c -lm -o boussinesqc; ./boussinesqc | gnuplot

note that using qcc it works like with cc.

Results

One analytical solution is the self similar solution \displaystyle \eta = \eta_0((t/2)^{-1/3}) Airy(2^{1/3} (x-t)t^{-1/3},t) the Airy function is approximated by its two asymptotic descirptions: \displaystyle sin(2/3.*((-x)^{3./2}+\pi/4))/\sqrt{(\pi*\sqrt{(-x)})}\text{ and }exp(-2/3 (x)^{(3./2)})/\sqrt{(4*\pi*\sqrt x)} and is in gnuplot as well airy

set xlabel "x"
plot [-20:]for [IDX=0:100:10] 'sol.OUT' i IDX u  1:2 not w lines

(script)

which gives U(x,t) plotted here for t=0 1 2 3 and -3<x<9

 dp3=2**(1./3)
 Airy(x) = (x<-.5? sin(2/3.*((-x)**(3./2)+pi/4))/sqrt(pi*sqrt(-x)): (x >.5 ? exp(-2/3.*(x)**(3./2))/sqrt(4*pi*sqrt(x)): NaN))
 set xlabel "(x-t)/t^{1/3}" 
 plot [-10:5]for [IDX=1:100:10] 'sol.OUT' i IDX u (($1-$3)/(($3/2)**(1./3))):($2*(($3/2)**(1./3))) not w lines, Airy(x),airy(x)