field.h

波浪数值模拟

/*
 * Copyright (c) 2001 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
 *
 */

/*  field.h

    C code to make and an array (field) data structure and do operations on it

*/


#ifndef FIELD_H
#define FIELD_H

#ifdef __cplusplus
extern "C" {
#endif				/* __cplusplus */


/* header file field.h  */

typedef enum {
  UGRID, VGRID, ETAGRID, UTGRID, VTGRID, VORTGRID, NOGRID    /* U=(nx,ny), V=(nx+1,ny), N=(nx-1,ny), UT=dU/dt=(nx-2,ny) VT=dV/dt=(nx-1,ny) */
} grid_t;

typedef struct  {
  int N;
  double *data;
} field1D;



typedef struct  {
  grid_t grid; 
  int N,M,num;
  double *data;
  double **row;
} field2D;

/* defines for dereferencing */

#define DR1(T,i)  ((T)->data[(i)])

#define DR2(T,i,j)  ((T)->data[(i)*((T)->M)+(j)])
#define DR2a(T,k)   ((T)->data[(k)])       // this is for when k=i*M+j

#define REF1(T)  ((T)->data)
#define REF2(T,i) ((T)->row[(i)])

  /* Defines for checking whether allocated */

#define NOT_ALLOCATED( T )  ( (T)==NULL )


/* field1D function declarations */

field1D*  allocate_field1D( int );
void deallocate_field1D(field1D * );
int  set_field1D_to_value(field1D *T, double a);

/* field 2D function declarations */

/* memory allocation & deallocation routines */

field2D*  allocate_field2D( int , int );
field2D*  allocate_field2D_grid( int , int , grid_t );
void deallocate_field2D(field2D * );
void deallocate_all_field2D();     /* use this before ending program to deallocate all field2D even the static ones */
double field2D_memory_usage();     /* Returns the number of megabytes allocated as field2D */

/* field2D arithmetic routinues:  First Addition then Multiplication */

int  field2D_set_to_value(field2D *T, double a);           /* T(:,:) = a  */
void field2D_copy(field2D *T, const field2D *A);           /* T = A */
void field2D_add(field2D *T, const field2D *A, const field2D *B);        /* T = A + B */
void field2D_add_sizes(field2D *T, const field2D *A, const field2D *B);  /* T = A + B : where A or B can be N+2 bigger than T */ 
void field2D_self_add(field2D *T, const field2D *A);       /* T += A  */
void field2D_self_add_two(field2D *T, const field2D *A, const field2D *B);     /* T += A+B  */
void field2D_self_add_three(field2D *T, const field2D *A, const field2D *B, const field2D *C);     /* T += A+B+C  */
void field2D_add_const(field2D *T, const double a);        /* T += ones(N,M)*a */
void field2D_mult_const(field2D *T, const double a);
void field2D_multiply(field2D *T, const field2D *A, const field2D *B);        /* T = A * B */
void field2D_multiply_sup1(field2D *T, const field2D *A, const field2D *B);   /* T = A * B */
void field2D_multiply_sub1(field2D *T, const field2D *A, const field2D *B);   /* T = A * B:  T(N,M), A(N-2,M), B(N,M) */
void field2D_self_multiply(field2D *T, const field2D *A);                     /* T *= A  */
void field2D_self_multiply_two(field2D *T, const field2D *A, const field2D *B);   /* T *= (A*B)  */
void field2D_multiply_add(field2D *T, const field2D *A, const field2D *B);    /* T  += A * B */

void field2D_self_divide(field2D *T, const field2D *A);                     /* T /= A  */

/* Input functions declarations */

int load_field1D_ascii(char *fname, field1D *T);    /* returns 1 if successfull */
int load_field1D_binary(char *fname, field1D *T);    /* returns 1 if successfull */

int load_field2D_ascii(char *fname, field2D *T);    /* returns 1 if successfull */
int load_field2D_binary(char *fname, field2D *T);    /* returns 1 if successfull */


/* Output function declarations */


void save_field1D_ascii(char *fname, field1D *T);
void save_field1D_binary(char *fname, field1D *T);
void print_field1D(field1D *T);

void print_field2D(field2D *T);
void save_field2D_ascii(char *fname, field2D *T);
void save_field2D_binary(char *fname, field2D *T);

void field_error_message(char *s);

  /* C-Grid operations for field2D */

  /* First the \bar operations */

void field2D_barx(field2D *FB, const field2D *F);
void field2D_bary_v_to_eta(field2D *ABY, const field2D *A);
void field2D_bary_eta_to_v(field2D *ABY, const field2D *A);
void field2D_bary_u_to_vort(field2D *ABY, const field2D *A);
void field2D_barxy_eta_to_vort(field2D *ABXY, const field2D *A);

  //void field2D_bary_vort_to_u(field2D *ABY, const field2D *A);

  /* Next the \delta_x and \delta_y operations AND gradient operator */

void field2D_diffx(field2D *DF, const field2D *F, double idx);
void field2D_diffx_add(field2D *DF, const field2D *F, double idx);
void field2D_diffy_eta_to_v(field2D *DF, const field2D *F, double idy);
void field2D_diffy_eta_to_v_add(field2D *DF, const field2D *F, double idy);
void field2D_diffy_v_to_eta(field2D *DF, const field2D *F, const double idy);
void field2D_diffy_v_to_eta_add(field2D *DF, const field2D *F, const double idy);
void field2D_diffy_u_to_vort(field2D *DF, const field2D *F, double idy);
void field2D_diffy_vort_to_u(field2D *DF, const field2D *F, double idy); 
void field2D_gradient_eta_to_uv(field2D *AX, field2D *AY, const field2D *A, double idx, double idy);

  /* Now 2nd derivative operations */

void field2D_diffx2(field2D *FXX, const field2D *F, const double idx);  /* \delta_x \delta_x (F) */
void field2D_diffy2(field2D *FYY, const field2D *F, const double idy);  /* \delta_x \delta_x (F) */
void field2D_laplacian(field2D *LAP_F, const field2D *F, const double idx, const double idy);  /* LAP_F = \delta_x^2 + \delta_y^2 (F) */ 

  /* Boundary condition operators */

void field2D_apply_nogradient_bc(field2D *T);     /* Apply dT/dx = 0 at boundary: i=0 and i=N */
void field2D_apply_zero_bc(field2D *T);           /* Apply T = 0 at boundary: i=0 and i=N     */

  /* generic debugging function */

double max_field2D(field2D * );
double min_field2D(field2D * );
double max_array(double *arr, int num);
int    lmax_array(double *arr, int num);

  /* Threaded functions */

  /* Initialization stuff */

void field2D_init_threads(int numcpu);
void field2D_close_threads();

  /* field2D_arithmetic_thread()  operations */

void field2D_copy_threads(field2D *T, const field2D *A);           /* T = A */
void field2D_self_add_threads(field2D *T, const field2D *A);       /* T += A  */
void field2D_add_const_threads(field2D *T, const double a);        /* T += ones(N,M)*a */
void field2D_mult_const_threads(field2D *T, const double a);
void field2D_multiply_threads(field2D *T, const field2D *A, const field2D *B);        /* T = A * B */
void field2D_multiply_sup1_threads(field2D *T, const field2D *A, const field2D *B);   /* T = A * B */
void field2D_self_multiply_threads(field2D *T, const field2D *A);                     /* T *= A  */
void field2D_multiply_add_threads(field2D *T, const field2D *A, const field2D *B);    /* T  += A * B */


  /* CGRID operations */


#ifdef __cplusplus
}                               /* extern "C" */
#endif                          /* __cplusplus */

#endif				/* FIELD_H */