field.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
 *
 */



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

#include <pthread.h>

#include "field.h"

typedef struct {
  int num_allocated;
  GSList *field2D_list;
} field2D_info_t;


field2D_info_t   field2D_info = {0, NULL };


void deallocate_all_field2D()
{
 
  GSList *list;

  for (list=field2D_info.field2D_list;list;list=list->next) 
    deallocate_field2D( (field2D *) list->data );

  g_slist_free(field2D_info.field2D_list);

}

/* Returns the number of megabytes allocated as field2D */
 
double field2D_memory_usage()
{

  GSList *list;
  field2D *T;
  int total=0;
  double megs;

  for (list=field2D_info.field2D_list;list;list=list->next) {
    T = (field2D * ) list->data;
    total += T->num;
  }

  //  fprintf(stderr,"
  megs = 8.0 * total/ (1024.0 * 1024.0);

  fprintf(stderr,"field2D:  number allocated =%d,  Memory requirements: %7.3f Megabytes\n",field2D_info.num_allocated,megs);

/* Writes out the estimated memory requirements for a model run 
   Assumes that almost all of the memory usage is due to field2D data structures */

void estimate_memory( )
{
    double mem;

    mem = field2D_memory_usage();  


}


  return megs;
}    

/* function definitions */


void field_error_message(char *s)
{

  fprintf(stderr,"** field.c: ERROR: %s\n",s);
  exit(-1);

}


/* field1D routinues */

field1D*  allocate_field1D(int N)
{
  field1D *T;

#ifdef __G_LIB_H_   /* if glib is defined */
  T = (field1D *) g_malloc(sizeof(field1D));
  T->data = (double *) g_malloc(N*sizeof(double));
#else
  T = (field1D *) malloc(sizeof(field1D));
  T->data = (double *) malloc(N*sizeof(double));
#endif  /* GLIB */

  T->N = N;

  return T;
}

field2D*  allocate_field2D(int N,int M)
{
  int i;
  field2D *T;

#ifdef __G_LIB_H_   /* if glib is defined */
  T = (field2D *) g_malloc(sizeof(field2D));  
  T->data = (double *) g_malloc(N*M*sizeof(double));
  T->row = (double **) g_malloc(N*sizeof(void*));
#else
  T = (field2D *) malloc(sizeof(field2D));  
  T->data = (double *) malloc(N*M*sizeof(double));
  T->row = (double **) malloc(N*sizeof(void*));
#endif

  T->grid = NOGRID;
  T->N = N;
  T->M = M;

  T->num = N*M;


  for (i=0;i<N;i++) 
    T->row[i] = &(T->data[i*M]);

  field2D_info.field2D_list = g_slist_append(field2D_info.field2D_list, T);
  field2D_info.num_allocated++;

  return T;
}


field2D*  allocate_field2D_grid(int N,int M, grid_t grid)
{
  int i;
  field2D *T;

#ifdef __G_LIB_H_   /* if glib is defined */
  T = (field2D *) g_malloc(sizeof(field2D));  
  T->data = (double *) g_malloc(N*M*sizeof(double));
  T->row = (double **) g_malloc(N*sizeof(void*));
#else  

  T = (field2D *) malloc(sizeof(field2D));  
  T->data = (double *) malloc(N*M*sizeof(double));
  T->row = (double **) malloc(N*sizeof(void*));
#endif

  T->grid = grid;
  T->N = N;
  T->M = M;

  T->num = N*M;


  for (i=0;i<N;i++) 
    T->row[i] = &(T->data[i*M]);


  /* Now add the field2D to the linked list for memory management purposes */

  field2D_info.field2D_list = g_slist_append(field2D_info.field2D_list, T);
  field2D_info.num_allocated++;

  /* eventually take this out right now just for testing */

  //  field2D_set_to_value(T, -8000.0);

  /* Return the pointer to the field2D */


  return T;
}

void deallocate_field1D(field1D *T)
{

  if (T!=NULL) {
    if (T->data!=NULL) {
#ifdef __G_LIB_H_   /* if glib is defined */
      g_free(T->data);
#else
      free(T->data);
#endif /* GLIB */
    }


#ifdef __G_LIB_H_   /* if glib is defined */
    g_free(T);
#else
    free(T);
#endif /* GLIB */
   }
}


void deallocate_field2D(field2D *T)
{

  if (T != NULL ) {
    if (T->row!=NULL)
#ifdef __G_LIB_H_   /* if glib is defined */
      g_free(T->row);
#else
      free(T->row);
#endif /* GLIB */


    if (T->data!=NULL) {
#ifdef __G_LIB_H_   /* if glib is defined */
      g_free(T->data);
#else
      free(T->data);
#endif /* GLIB */
    }

#ifdef __G_LIB_H_   /* if glib is defined */
    g_free(T);
#else
    free(T);
#endif /* GLIB */
   }

}



int  set_field1D_to_value(field1D *T, double a)
{
  int i;

  if ((T==NULL)||(T->data==NULL))
    return 0;

  for (i=0;i<T->N;i++)
    T->data[i] = a;

  return 1;
}

int  field2D_set_to_value(field2D *T, double a)
{
  int i;

  double *td;

  if ((T==NULL)||(T->data==NULL))
    return 0;

  td = T->data;

  for (i=0;i<T->num;i++)
    *td++  = a;

  return 1;
}


void field2D_mult_const(field2D *T, const double a)
{
  int i;
  double *td = T->data;

  for (i=0;i<T->num;i++)
    *td++ *= a;

}


double max_field2D(field2D *T)
{
  int i;
  double tmp,max=-1.0;

  for (i=0;i<T->num;i++) {
    tmp = fabs(T->data[i]);
    if (tmp>max) 
      max = tmp;
  }

  return max;
}


double min_field2D(field2D *T)
{
  int i;
  double tmp,min=1e+38;

  for (i=0;i<T->num;i++) {
    tmp = fabs(T->data[i]);
    if (tmp<min) 
      min = tmp;
  }

  return min;
}


/* Field2D math routinues */

/* Copy routine, T = A:  Dimensions:  T(N,M) and A is either A(N,M) and A(N-2,M)  */

void field2D_copy(field2D *T, const field2D *A)               /* T = A */
{
  int NUM = T->num;
  int N = T->N;
  int M = T->M;
  int i;

  double *td = T->data;
  const double *ad = A->data;

  /* Check that the dimensions are ok */

  if ( A->M != M  ) {
    field_error_message("Bad field2D M sizes to void field2D_copy()");
  }

  if ( !( ( A->N == N) || (A->N == N-2) ) ) {
    field_error_message("Bad field2D N sizes to void field2D_copy()");
  }

  switch (A->N - N) {
  case 0:           /* if A & T have same sizes, OK, start at same location */
    ad = A->data;
    td = T->data;
    NUM = T->num;
    break;
  case -2:             /* if A(N,M) & T(N+2,M) then move td pointer up */
    ad = A->data;
    td = &(T->data[M]);
    NUM = A->num;
    break;
  default:  field_error_message("Bad case to switch() {...}  in field2D_copy()");
  }

  /* Then go through and do the copy */

  for (i=0;i<NUM;i++)
    *td++ = *ad++;

}


/* Field2D multiplly routines:   T = A + B */

void field2D_add(field2D *T, const field2D *A, const field2D *B)
{

  int NUM = T->num;
  int N = T->N;
  int M = T->M;
  int i;

  double *td = T->data;
  const double *ad = A->data;
  const double *bd = B->data;

  /* Check that the dimensions are ok with an #ifdef */

  if ( (A->N != N) || (B->N != N) || (A->M != M) || (B->M != M) ) {
    field_error_message("Bad field2D sizes to void field2D_add()");
  }

  /* Then go through and do the multiplication */

  for (i=0;i<NUM;i++)
    *td++ = (*ad++) + (*bd++);

}


/* Field2D multiplly routines:   T = A + B   
    - this version allows for mismatched sizes.   So
    either  A or B (or both) can be n+2 bigger than T.   */

void field2D_add_sizes(field2D *T, const field2D *A, const field2D *B)
{

  int NUM = T->num;
  int N = T->N;
  int M = T->M;
  int i;

  double *td = T->data;
  const double *ad = A->data;
  const double *bd = B->data;

  /* Check that the dimensions are ok with an #ifdef */

  if ( (A->M != M) || (B->M != M) ) {
    field_error_message("Bad M field2D sizes to void field2D_add_sizes()");
  }



  switch( A->N - N ) {
  case 0: ad = A->data;
    break;
  case 2: ad = REF2(A, 1);
    break;
  default: field_error_message("Bad N field2D sizes for A in field2D_add_sizes()");
  }

  switch( B->N - N ) {
  case 0: bd = B->data;
    break;
  case 2: bd = REF2(B, 1);
    break;
  default: field_error_message("Bad N field2D sizes for B in field2D_add_sizes()");
  }

  /* Then go through and do the multiplication */

  for (i=0;i<NUM;i++)
    *td++ = (*ad++) + (*bd++);

}



/* Field2D multiplly routines:   T = A * B */

void field2D_multiply(field2D *T, const field2D *A, const field2D *B)
{

  int NUM = T->num;
  int N = T->N;
  int M = T->M;
  int i;

  double *td = T->data;
  const double *ad = A->data;
  const double *bd = B->data;

  /* Check that the dimensions are ok with an #ifdef */

  if ( (A->N != N) || (B->N != N) || (A->M != M) || (B->M != M) ) {
    field_error_message("Bad field2D sizes to void field2D_multiply()");
  }

  /* Then go through and do the multiplication */

  for (i=0;i<NUM;i++)
    *td++ = (*ad++) * (*bd++);

}

/* Field2D multiplly routines:   T = A * B */
/* Assumes  T = (N,M),   A = (N+2,M),  B = (N,M)  - for A multiplies in the middle  A(2:N+1,:) */

void field2D_multiply_sup1(field2D *T, const field2D *A, const field2D *B)
{

  int NUM = T->num;
  int N = T->N;
  int M = T->M;
  int i;

  double *td = T->data;
  const double *ad = &(A->data[M]);
  const double *bd = B->data;

  /* Check that the dimensions are ok with an #ifdef */

  if ( (A->N != N+2) || (B->N != N) || (A->M != M) || (B->M != M) ) {
    field_error_message("Bad field2D sizes to void field2D_multiply()");
  }

  /* Then go through and do the multiplication */

  for (i=0;i<NUM;i++)
    *td++ = (*ad++) * (*bd++);

}


/* Field2D multiplly routines:   T = A * B */
/* Assumes  T = (N,M),   A = (N-2,M),  B = (N,M)  - for A multiplies in the middle  A(2:N+1,:) */

void field2D_multiply_sub1(field2D *T, const field2D *A, const field2D *B)
{


  int N = T->N;
  int M = T->M;
  int NUM = T->num - 2*M;

  int i;

  double *td = REF2(T, 1);
  const double *ad = REF2(A, 0);
  const double *bd = REF2(B, 1);

  /* Check that the dimensions are ok with an #ifdef */

  if ( (A->N != N-2) || (B->N != N) || (A->M != M) || (B->M != M) ) {
    field_error_message("Bad field2D sizes to void field2D_multiply_sub1()");
  }

  /* Then go through and do the multiplication */

  for (i=0;i<NUM;i++)
    *td++ = (*ad++) * (*bd++);

}


/* Field2D math routinues :   T += A  
    If T has size (N,M), then A must also be size (N,M) OR
    - if A is (N+2,M) then the multiply is done from (2:N-1,M)  */

void field2D_self_add(field2D *T, const field2D *A)
{

  int N = T->N;
  int M = T->M;
  int i;
  int NUM;

  double *td = T->data;
  const double *ad = A->data;

  /* Check that the dimensions are ok with an #ifdef */

  if ( (A->M != M) ) { 
    field_error_message("Bad M field2D sizes to void field2D_self_add()");
  }

  switch (A->N - N) {
  case 0:       /* if A & T have same sizes, OK, start at same location */
    ad = A->data;
    td = T->data;
    NUM = T->num;
    break;
  case 2: 
    ad = &(A->data[M]); /* if A(N+2,M) & T(N,M) then move ad pointer up */
    td = T->data;
    NUM = T->num;
    break;
  case -2:             /* if A(N,M) & T(N+2,M) then move td pointer up */
    ad = A->data;
    td = &(T->data[M]);
    NUM = A->num;
    break;
  default:   field_error_message("Bad N field2D sizes to void field2D_self_add()");
  }


  /* Then go through and do the addition */

  for (i=0;i<NUM;i++) {
    *td++ += *ad++; 
  }
}


/* Field2D math routinues :   T += A+B+C