floats.c

波浪数值模拟

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



/* floats.c    Falk Feddersen */       

#include <stdlib.h>      /* this is for the function rand() */
#include <time.h>        /* this is to get the seed for random number generator */
#include <math.h>        /* need this for the ceil() and floor() functions */

#include "floats.h"

floats_t  Floats;

switch_t  floats_flag;

static FILE *outfp;

double Ly, Lxlimit, L0limit;

void free_floats()
{

  if (outfp!=NULL) {   /* If the outfile isn't closed then close it here */
    fclose(outfp);
  }

  g_free(Floats.xloc);
  g_free(Floats.yloc);
  g_free(Floats.uvel);
  g_free(Floats.vvel);

}

void floats_report(floats_info_t *FL)
{

  fprintf(stdout,"------------- * Floats Report * --------------------------\n");

  /* Then model dimension information */
  /* - Passive Tracer - */

  if (FL->floats_flag == ON) {
    switch (FL->floatinit) {
    case RANDOM_FLOATS:  fprintf(stdout,"Floats: Random Floats, num=%d\n",FL->nfloats);
      break;
    case FILE_FLOATS:  fprintf(stdout,"Floats: File Floats, file=%s,  num=%d\n",FL->floats_file->str,Floats.nfloats);
      break;
    case POINT_FLOATS: fprintf(stdout,"Floats: Point Floats, num=%d\n",Floats.nfloats);
      break;
    }
  }
  else
    fprintf(stdout,"Floats:\t module not used in this run\n");

  fprintf(stdout,"------------------------------------------------------------------\n");
}


void init_random_floats(floats_info_t *FL)
{

  time_t tt;
  unsigned int seed;
  int i;

  Floats.nfloats = FL->nfloats;
  if (Floats.nfloats>0) {
    Floats.xloc = (double * ) g_malloc(sizeof(double)*Floats.nfloats);
    Floats.yloc = (double * ) g_malloc(sizeof(double)*Floats.nfloats);
    Floats.uvel = (double * ) g_malloc(sizeof(double)*Floats.nfloats);
    Floats.vvel = (double * ) g_malloc(sizeof(double)*Floats.nfloats);
  }

  /* seed the random number generator */

  time(&tt);
  seed = (int ) tt;
  srand(seed);

  /* generate the random values */

  for (i=0;i<Floats.nfloats;i++) {
    Floats.xloc[i] =  nx*dx* (((double ) rand())/RAND_MAX);
    Floats.yloc[i] =  ny*dy* (((double ) rand())/RAND_MAX);
  }
}


/* The floats input file format is rows of  x & y  (two columns)
   The strategy for inputing is to read the file once to get the total number
   of floats.   Then to create the dynamic arrays and then re-read the floasts file
   to get the x & y positions */

void load_floats_info(GString *float_file)
{
  int num, i;
  char s[120];
  double xx, yy;
  FILE *fp;

  if ((fp = fopen(float_file->str,"r"))==NULL) {
    funwaveC_error_message("Cannot open input floats file");
  }

  /* Read the number of floats */

  while (fgets(s,120,fp)!=NULL) {
    num++;
  }

  fclose(fp);

  /* save the number of floats */
  Floats.nfloats = num;

  /* re-open the floats file */

  if ((fp = fopen(float_file->str,"r"))==NULL) {
    funwaveC_error_message("Cannot open input floats file");
  }

  for (i=0;i<Floats.nfloats;i++) {
    fscanf(fp,"%lf %lf",&xx, &yy);
    Floats.xloc[i] = xx;
    Floats.yloc[i] = yy;
  }

  fclose(fp);



}


void init_floats(floats_info_t *floatsinfo)
{
  int num;

  num = floatsinfo->nfloats;
  
  floats_flag = floatsinfo->floats_flag;

  Lxlimit = (nx-1)*dx - 0.25*dx;
  L0limit = 0.25*dx;
  Ly = ny*dy;

  if (num<0)
    funwaveC_error_message("Less than zero floats specified!");

  switch (floatsinfo->floatinit ) {
  case FILE_FLOATS:  load_floats_info(floatsinfo->floats_file);
    break;
  case RANDOM_FLOATS:  init_random_floats(floatsinfo);
    break;
  }


}




/* Function to calculate the float velocities by interpolation over the grid */

void calculate_float_velocities(const field2D *UU, const field2D *VV)
{

  int k;

  int ix1, ix2, jx1, jx2;
  int iy1, iy2, jy1, jy2;

  double xx, yy, xx2, yy2;
  double px1, px2, py1, py2;

  double *uud = Floats.uvel;
  double *vvd = Floats.vvel;

  double dxdy = dx/dy;
  double idxdy = idx/idy;

  for (k=0;k<Floats.nfloats;k++) {

    xx = Floats.xloc[k]*idx;  /* xx & yy are float positions normalized by the grid size */
    yy = Floats.yloc[k]*idy;   

    /* First do the U velocities */

    ix1 = (int ) floor(xx);
    ix2 = (int ) ceil(xx);   // so ix2 should be ix1 + 1 (but may not be)
    jx1 = (int ) floor(yy);
    jx2 = (int ) ceil(yy);

    if (jx2>(ny-1))   /* allow this to wrap around for velocity interpolation */
      jx2 = 0;    

    px2 = xx-ix1;
    px1 = 1.0-px2;
    py2 = yy-jx1;
    py1 = 1.0 - py2;

    *uud++ =  px1 * ( py1*DR2(UU,ix1,jx1) + py2*DR2(UU,ix1,jx2) ) +  
              px2 * ( py1*DR2(UU,ix2,jx1) + py2*DR2(UU,ix2,jx2) );


    /* Next do the V velocities */

    xx2 = xx+0.5;
    yy2 = yy-0.5;
    iy1 = (int ) floor(xx2);
    iy2 = (int ) ceil(xx2);   // so ix2 should be ix1 + 1 (but may not be)
    jy1 = (int ) floor(yy2);
    jy2 = (int ) ceil(yy2);

    if (jy1<0)      /* take care of the possibilities that we are at the boundary */ 
      jy1 = ny-1;

    if (jy2>(ny-1)) 
      jy2 = 0;
    
    px2 = xx2-iy1;
    px1 = 1.0-px2;
    py2 = yy2-jy1;
    py1 = 1.0 - py2;

    *vvd++ =  px1 * ( py1*DR2(VV,iy1,jy1) + py2*DR2(VV,iy1,jy2) ) +  
              px2 * ( py1*DR2(VV,iy2,jy1) + py2*DR2(VV,iy2,jy2) );

  }
              

}


/* Function that uses a first order time stepping to move the floats */


void step_floats_one(double dt)   
{

  int i;
  double xnew, ynew;

  double *uud = Floats.uvel;
  double *vvd = Floats.vvel;
  double *xd = Floats.xloc;
  double *yd = Floats.yloc;
  int  N = Floats.nfloats;

  for (i=0;i<N;i++) {
    xnew = *xd + (*uud++  * dt);
    ynew = *yd + (*vvd++  * dt);

    if (ynew > Ly)
      ynew -= Ly;
    if (ynew <0)
      ynew +=  Ly;

    *yd++ = ynew;

    if (xnew > Lxlimit)
      xnew = Lxlimit;
    if (xnew < L0limit)
      xnew = L0limit;

    *xd++ = xnew;

  }


}


/* 
   The algorithm for this is this:
     => for each float
          - calculate the nearest 4 U & V locations 
          - interpolate them onto the float location
          - update w/ a single euler step the float velocity

      ? How does ROMS do it ?
*/

   

void update_floats_AB2(const field2D *UU, const field2D *VV, double dt)
{

  int i;

  calculate_float_velocities(UU, VV);
  step_floats_one(dt);

}


void update_floats_euler(const field2D *UU, const field2D *VV, double dt)
{

  int i;

  calculate_float_velocities(UU, VV);
  step_floats_one(dt);

}