tridiagonal.c

波浪数值模拟

/*
 * Copyright (c) 2001-2005 Falk Feddersen
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */



/* tridiagonal.c   Falk Feddersen */       

#include <assert.h>
#include "tridiagonal.h"
//#include "bousinessq_constants.h"


/* /\* #ifdef LINEAR *\/ */
/* /\* const double trib1 = 0.0; //0.5* beta*beta;        /\\* trib1 and b2 are bousinessq parameters in the d(u,v)/dt = ... equation *\\/ *\/ */
/* /\* const double trib2 = 0.0; //beta; *\/ */
/* /\* #endif *\/ */

/* /\* #ifdef NSWE *\/ */
/* /\* const double trib1 = 0.0; //0.5* beta*beta;        /\\* b1 and b2 are bousinessq parameters in the d(u,v)/dt = ... equation *\\/ *\/ */
/* /\* const double trib2 = 0.0; //beta; *\/ */
/* /\* #endif *\/ */

/* /\* #ifdef PEREGRINE   /\\* peregrine depth-averaged dynamics *\\/ *\/ */
/* /\* const double trib1 = 1.0/6.0; /\\* not sure why b1 is what it is *\\/ *\/ */
/* /\* const double trib2 = -0.5;    /\\* or b2 but these are the peregrine JFM (1967) values *\\/ *\/ */
/* /\* #endif  /\\* PEREGRINE *\\/ *\/ */

/* /\* #ifdef WEI_KIRBY *\/ */
/* /\* const double trib1 =   1.409805000000000e-01;    /\\* b1 = 0.5* beta_h*beta_h;   b1 and b2 are bousinessq parameters in the d(u,v)/dt = *\\/ *\/ */
/* /\* const double trib2 =  -0.531; *\/ */
/* /\* #endif *\/ */

/* /\* #ifdef NWOGU *\/ */
/* /\* const double trib1 =   1.409805000000000e-01;    /\\* b1 = 0.5* beta_h*beta_h;   b1 and b2 are bousinessq parameters in the d(u,v)/dt = ... equation *\\/ *\/ */
/* /\* const double trib2 =  -0.531; *\/ */
/* /\* #endif *\/ */

double trib1;
double trib2;

/* Just make these global variables.  No other .c or .h file sees them */

/* The U tridiagonal solve variables */

field2D *bU;     /* the transpose of RHS_U */
field2D *UTMP;   /* the result of tridiagonal_solve_u - transpose to get U */
field2D *AD_U;   /* The diagonal at each j for U*/
field2D *ADL_U;  /* The lower diagonal at each j for U */
field2D *ADU_U;  /* the upper diagonal at each j for U */

field2D *AD_V;   /* The diagonal at each j for V */
field2D *ADL_V;  /* The lower diagonal at each j for V  ADL_V(i,1) is the A(1,n) term in the tridiagonal matrix */
field2D *ADL2_V;  /* using the GNU scientific library routinue */
field2D *ADU_V;  /* the upper diagonal at each j for V  ADU_V(i,MV-1) is the A(n,1) term in the tridiagonal matrix */


//double *b1v;     /* The b(1) term in the V periodic tri-diag matrix - at every i index */
//double *bnv;     /* The b(n) term in the V periodic tri-diag matrix - at every i index */

double *worku;
double *workv;  /* work arrays for tridiagonal solving */

void  finish_init_tridiagonal_v(int NV, int MV );

void  make_ADL2_from_ADL(field2D *ADL2_V, const field2D *ADL);


int check_v_tridiagonal_solution(double *veq_ad, double *veq_adl, double *veq_adu, double *vd, double *rhs_vd, int M);
int check_v_vtmp_integral(double *vd, double *rhs_vd, int M, int index);



void  tri_init_bousinessq_constants(dynamics_t Dynamics)
{
  switch (Dynamics) {
  case WEI_KIRBY_DYNAMICS:
    trib1 =   1.409805000000000e-01;    /* trib1 = 0.5* beta_h*beta_h;   trib1 and trib2 are bousinessq parameters in the d(u,v)/dt = */
    trib2 =  -0.531;
    break;
  case NWOGU_DYNAMICS:
    trib1 =   1.409805000000000e-01;    /* trib1 = 0.5* beta_h*beta_h;   trib1 and trib2 are bousinessq parameters in the d(u,v)/dt = ... equation */
    trib2 =  -0.531;
    break;
  case PEREGRINE_DYNAMICS:
    trib1 = 1.0/6.0; /* not sure why trib1 is what it is */
    trib2 = -0.5;    /* or trib2 but these are the peregrine JFM (1967) values */
    break;
  case NSWE_DYNAMICS:
    trib1 = 0.0; //0.5* beta*beta;        /* trib1 and trib2 are bousinessq parameters in the d(u,v)/dt = ... equation */
    trib2 = 0.0; //beta;
    break;
  case LINEAR_DYNAMICS:
    trib1 = 0.0; //0.5* beta*beta;        /* trib1 and trib2 are bousinessq parameters in the d(u,v)/dt = ... equation */
    trib2 = 0.0; //beta;
    break;
  }
}

void  init_tridiagonal_solver(field2D *HBX, field2D *HBY, dynamics_t Dynamics)
{
  int NV = HBY->N;
  //  int NV2 = NV-2;
  int MV = HBY->M;
  int MV2 = MV-2;

  int NU = HBX->N;   /* x size of U */
  int NU2 = NU-2;    /* NU - 2 for ADU and ADL and AD  */
  int MU = HBX->M;


  int i,j;
  double hbxi, hbxip, hbxim;
  double hbyi, hbyip, hbyim;
  double tmpa0, tmpb0, tmp1d;      /* for the main diagonal */

  /* First initialize the constants */

  tri_init_bousinessq_constants(Dynamics);

  /* allocate the U fields - ADL & ADU fields etc. - not sizes are transposed! */

  UTMP = allocate_field2D(MU,NU);
  bU = allocate_field2D(MU,NU);    /* this is the transposed rhs of U eq */
  AD_U = allocate_field2D(MU,NU2);
  ADL_U = allocate_field2D(MU,NU2);
  ADU_U = allocate_field2D(MU,NU2);

  /* allocate the V memeory - first the V ADL & ADU fields etc.  */

  AD_V = allocate_field2D(NV,MV);
  ADL_V = allocate_field2D(NV,MV);
  ADL2_V = allocate_field2D(NV,MV);
  ADU_V = allocate_field2D(NV,MV);

  /* the work spaces are next */

  worku = (double * ) g_malloc(NU2*sizeof(double));
  workv = (double * ) g_malloc(MV2*sizeof(double));

  /* Now initialize these guys with the HBX and HBY inputs */
  /* first do U */

  for (j=0;j<MU;j++) {
    for (i=0;i<NU2;i++) {
      hbxi = DR2(HBX,i+1,j);
      hbxip = DR2(HBX,i+2,j);
      hbxim = DR2(HBX,i,j);
      tmpa0 = trib1*idx2*hbxi*hbxi;
      tmpb0 = trib2*idx2*hbxi*hbxi;
      tmp1d = -2.0*(tmpa0 + tmpb0);
      DR2(AD_U,j,i) = 1 + tmp1d;
      DR2(ADU_U,j,i) = tmpa0 + trib2*idx2*hbxi*hbxip;
      DR2(ADL_U,j,i) = tmpa0 + trib2*idx2*hbxi*hbxim;
    }
    DR2(ADL_U,j,0) = 0.0;    /* make these last guys be zeros!! */
    DR2(ADU_U,j,NU2-1) = 0.0;
  }

  /* This sets up the periodic tridiagonal system for V - the upper, lower, and main diagonals */
  /* Need to solve  v + h [b1*h*v_yy + b2 (hv)_yy ] = rhs
  /* only loop through the interior i points because periodic V bc will be applied */

  for (i=1;i<NV-1;i++) { 

    /* First do j = 0 */
    j = 0;

    hbyi = DR2(HBY,i,j);       
    hbyip = DR2(HBY,i,j+1);
    hbyim = DR2(HBY,i,MV-1);
    tmpa0 = trib1*idy2*hbyi*hbyi;
    tmpb0 = trib2*idy2*hbyi*hbyi;
    tmp1d = -2.0*(tmpa0 + tmpb0);
    DR2(AD_V,i,j) = 1 + tmp1d;
    DR2(ADU_V,i,j) = tmpa0 + trib2*idy2*hbyi*hbyip;
    DR2(ADL_V,i,j) = tmpa0 + trib2*idy2*hbyi*hbyim;  // ADL_V(i,0) is actually the A(1,n) term in the tridiagonal matrix

    for (j=1;j<MV-1;j++) {
      hbyi = DR2(HBY,i,j);
      hbyip = DR2(HBY,i,j+1);
      hbyim = DR2(HBY,i,j-1);
      tmpa0 = trib1*idy2*hbyi*hbyi;
      tmpb0 = trib2*idy2*hbyi*hbyi;
      tmp1d = -2.0*(tmpa0 + tmpb0);
      DR2(AD_V,i,j) = 1 + tmp1d;
      DR2(ADU_V,i,j) = tmpa0 + trib2*idy2*hbyi*hbyip;
      DR2(ADL_V,i,j) = tmpa0 + trib2*idy2*hbyi*hbyim;
    }

    /* Last do j = MV-1 */

    j = MV-1;

    hbyi = DR2(HBY,i,j);
    hbyip = DR2(HBY,i,0);
    hbyim = DR2(HBY,i,j-1);
    tmpa0 = trib1*idy2*hbyi*hbyi;
    tmpb0 = trib2*idy2*hbyi*hbyi;
    tmp1d = -2.0*(tmpa0 + tmpb0);
    DR2(AD_V,i,j) = 1 + tmp1d;
    DR2(ADU_V,i,j) = tmpa0 + trib2*idy2*hbyi*hbyip;  // ADL_V(i,MV-1) is actually the A(n,1) term in the tridiagonal matrix
    DR2(ADL_V,i,j) = tmpa0 + trib2*idy2*hbyi*hbyim;  // ADL_V(i,0) is actually the A(1,n) term in the tridiagonal matrix

  }

  make_ADL2_from_ADL(ADL2_V, ADL_V);

  finish_init_tridiagonal_v(NV, MV );
}

void  deallocate_tridiagonal_solver()
{

  deallocate_field2D(bU);
  deallocate_field2D(UTMP);
  deallocate_field2D(AD_U);
  deallocate_field2D(ADL_U);
  deallocate_field2D(ADU_U);

  deallocate_field2D(AD_V);
  deallocate_field2D(ADL_V);
  deallocate_field2D(ADU_V);

  g_free(worku);
  g_free(workv);

}

/* ------ TRIDIAGONAL MATRIX SOLVER ROUTINES ------------- 

tridiagonal matrix solver routinues
takes the diagonal, and the upper and lower and returns the answer in x
AD, ADU, and ADL are all length N vectors.

The first element of ADL is zero and the last element of ADU is zero.
on return:  0:  OK,     1:  BAD  (AD=0.0),   2:  ZERO PIVOT           
 !!! I could speed this up by pre setting up the tridiagonal solution 
  */

void  tridiagonal_solve_u(double *ad,double *adl,double *adu,double *bu,double *u,double *work,int n)
{

    int j;
    double bet;
    double *x = &(u[1]);
    double *b = &(bu[1]);

    assert(  (ad[0]!=0.0) );

    
    x[0]=b[0]/(bet=ad[0]);

    for (j=1;j<n;j++) {
       work[j]=adu[j-1]/bet;
       bet = ad[j]-adl[j]*work[j];
       if (bet == 0.0) {
	 funwaveC_error_message("** Bad result (bet==0.0) in tridiagonal.c: tridiagonal_solve_u()");
       }

       x[j]=(b[j]-adl[j]*x[j-1])/bet;
    }
    for (j=n-2;j>=0;j--) 
       x[j] -= work[j+1]*x[j+1];
    
    u[0] = 0.0;    /* Here the no normal flow conditions are applied */
    u[n+1] = 0.0;


}

/*  General solver for the V periodic tridiagonal solver that gets called repeatedly */

void  tridiagonal_solve_for_v_work(double *ad,double *adl,double *adu,double *b,double *x,double *work,int n)
{

    int j;
    double bet;

    assert( (ad[0]!=0.0) );
    
    x[0]=b[0]/(bet=ad[0]);

    for (j=1;j<n;j++) {
       work[j]=adu[j-1]/bet;
       bet = ad[j]-adl[j]*work[j];
       if (bet == 0.0) {
	 funwaveC_error_message_file_line("** Bad result (bet==0.0) in tridiagonal_solve_for_v_work()",__FILE__, __LINE__);
       }

       x[j]=(b[j]-adl[j]*x[j-1])/bet;
    }
    for (j=n-2;j>=0;j--) 
       x[j] -= work[j+1]*x[j+1];
    
}

// THIS is taken from periodic_tridiag.m
//
// The idea here is to write the matrix problem as follows:
//
// a11 v1 a1n   | x1       b1 
//              |
//  v2 T  v3    | xx    =  bb 
//              |
// an1 v4 ann   | xn       bn
//
// where v1,v2,v3,v4,xx,and bb are vectors of length n-2.  and T is a
// tridiagonal matrix.
//
// The problem then becomes
//
// a11 x1 + <v1,xx> + a1n xn = b1
// x1 v2 + T xx + xn v3 = bb
// an1 x1 + <v4,xx> + ann xn = bn
//
// We solve the second equation for xx:  
// T xx = bb - x1 v2 - xn v3
// xx = inv(T)*bb - x1 inv(T)*v2 - xn inv(T)*v3
//
// plug xx into the first and third equations, resulting in a 2x2 problem
// for x1 and xn
// 
// a11 x1 + <v1,inv(T)*bb> -x1 <v1,inv(T)*v2> - xn <v1,inv(T)*v3> + a1n xn = b1
// an1 x1 + <v4,inv(T)*bb> -x1 <v4,inv(T)*v2> - xn <v4,inv(T)*v3> + ann xn = bn
//
// AA(1,1) = a11 - <v1,inv(T)*v2>
// AA(1,2) = a1n - <v1,inv(T)*v3>
// AA(2,1) = an1 - <v4,inv(T)*v2>
// AA(2,2) = ann - <v4,inv(T)*v3>
//
// rhs(1) = b1-<v1,inv(T)*bb>
// rhs(2) = bn-<v4,inv(T)*bb>
//
// once we know x1 and xn, we use them to define xx:
// xx = inv(T)*bb - x1 inv(T)*v2 - xn inv(T)*v3
//
// Note: this uses three calls to tridiag, to find inv(T)*bb,inv(T)*v2,inv(T)*v3
// which is 3*O(n).  But this still beats the LU approach which is O(n^2).

field2D *V2i;     /* at each i or x that is tridiagonally solved, this is v2i(1:M-2) */
field2D *V3i;    /* at each i or x that is tridiagonally solved, this is v2i(1:M-2) */


void  finish_init_tridiagonal_v(int NV, int MV )
{
  int i, j;
  double *v2, *v3;
  double *tmp_array;
  double *veq_ad,  *veq_adl, *veq_adu;

  /* First allocate the temporary memory */

  v2 = (double * ) g_malloc((MV-2)*sizeof(double));
  v3 = (double * ) g_malloc((MV-2)*sizeof(double));
  tmp_array = (double * ) g_malloc((MV-2)*sizeof(double));
  for (i=0;i<MV-2;i++) {
    v2[i] = 0.0;
    v3[i] = 0.0;
    tmp_array[i] = 0.0;
  }

  V2i = allocate_field2D(NV-2,MV-2);
  V3i = allocate_field2D(NV-2,MV-2);


  /* Then solve each tridiagonal case at the various i 
     Only need to save v2i and v3i at the start and end of the array
     because it gets multiplied by the 1st and last element of v1 and v4 
     
     Note that only first and last elements of v2 and v3 are non-zero    */
 
  for (i=1;i<NV-1;i++) {
    v2[0] = DR2(ADL_V, i, 1);
    v3[MV-3] = DR2(ADL_V, i, MV-2);
    veq_ad = REF2(AD_V, i ) ; //&(AD_V->data[i*MVV]);
    veq_adl = REF2(ADL_V, i ); //&(ADL_V->data[i*MVV]);
    veq_adu = REF2(ADU_V, i ); // &(ADU_V->data[i*MVV]);

    tridiagonal_solve_for_v_work(veq_ad+1, veq_adl+1, veq_adu+1, v2, tmp_array , workv, MV-2);  

    for (j=0;j<MV-2;j++) {
      DR2(V2i,i-1,j) = tmp_array[j];
    }

    tridiagonal_solve_for_v_work(veq_ad+1, veq_adl+1, veq_adu+1, v3, tmp_array , workv, MV-2);  
    for (j=0;j<MV-2;j++) {
      DR2(V3i,i-1,j) = tmp_array[j];
    }

  }

  g_free(v2);
  g_free(v3);
}