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



/* initial_condition.c   Falk Feddersen */       

#include <stdlib.h>

#include "initial_condition.h"

/* This function should be called in the initialization phase to zero out
   all the matrices  */


/* Function to load the Initial Condition file.
   The flag is set to 1 for U, 2 for V

   We can adjust if we want IC's for P too */

void load_full_IC(char *fname, field2D *TT)
{

  int val;
  val = load_field2D_ascii(fname,TT);
  if (val != 1)
    funwaveC_error_message("Error loading the full IC file!\n");

}



/* load the homogeneous intitial condition for U,V, or P 

   */


void load_line_IC(char *fname, field2D *TT)
{

   FILE *fp;
   int i,j,up,ac;
   gdouble a;

   if ((fp=fopen(fname,"r"))==NULL) {
       fprintf(stderr,"Error loading line_IC file: %s\n",fname);
       exit(-1);
   }

   up = TT->N;
   ac = TT->M;

   for (i=0;i<up;i++) {
      fscanf(fp,"%lf",&a);
      for (j=0;j<ac;j++)
        DR2(TT,i,j) = a;
   }
   fclose(fp);

}





/* This function creates the y-forcing to balance
   the V initial condition before the perturbation
   is applied  */
/*
void balance_vic(field2D *FY, field2D *VV)
{
   int i,j;

   for (i=1;i<(VV->N)-1;i++)
     for (j=0;j<FY->M;j++)
        DR2(FY,i-1,j) = mu*DR2(VV,i,j);
  

}
*/



/* This sets the V IC to a constant velocity */


void set_const_v(field2D *VV, gdouble a)
{
#ifdef NOSLIP
  int j;
#endif  /* NOSLIP */

  field2D_set_to_value(VV,a);

#ifdef NOSLIP
   for (j=0;j<VV->M;j++) {
        DR2(VV,0,j) = -DR2(VV,1,j);
        DR2(VV,(VV->N)-1,j) = -DR2(VV,(VV->N)-2,j);
   }
#endif

}



/* this function reads in a perturbation streamfunction file and sets up the
   perturbation u & v velocities.    */

/* THIS FUNCTION IS BROKEN!!!! */
void  load_streamfunction(char *pname, field2D *U, field2D *V, field2D *H)
{
  int i,j;
  gdouble pp,t0,t1,hbar, depth;
  field2D *SF;

  int nx = U->N;
  int ny = U->M;

  fprintf(stderr,"** WARNING: in broken function input.c;load_streamfunction()!!");
  depth=1.0;

  SF = allocate_field2D(nx,ny);

  if (!load_field2D_ascii(pname,SF)) 
     funwaveC_error_message("Couldn't load streamfunction as field2D\n");

  /* Test the streamfunction to make sure that SF[0][j] & SF[NX-1][j] are homogeneous */

   t0 = DR2(SF,0,0);
   t1 = DR2(SF,nx-1,0);
   for (j=1;j<ny;j++) 
     if ((DR2(SF,0,j)!=t0)||(DR2(SF,nx-1,j)!=t1)) 
           funwaveC_error_message("***** Error in field2D streamfunction at boundary!  EXIT!");
 
   /* now generate the u and v perturbations */
   /* first do the V perturbations */

   for (j=0;j<ny;j++) {
      pp = (DR2(SF,1,j)-DR2(SF,0,j))*idx/depth;
      DR2(V,1,j) += pp;
      DR2(V,0,j) = DR2(V,1,j);

      for (i=1;i<nx-1;i++) {
         pp = (DR2(SF,i+1,j)-DR2(SF,i,j))*idx/depth;
         DR2(V,i+1,j) += pp;
      }
      DR2(V,nx,j) = DR2(V,nx-1,j);      

#ifdef NOSLIP
      DR2(V,0,j) = -DR2(V,1,j);
      DR2(V,nx,j) = -DR2(V,nx-1,j);      
#endif
   }

   /* next do the U perturbations */

   for (i=1;i<nx-1;i++) {
     hbar = -1.0;  //HB[i-1];
      pp = -(DR2(SF,i,0)-DR2(SF,i,ny-1))*idy/hbar;
      DR2(U,i,0) += pp;
      for (j=1;j<ny;j++) {
         pp = -(DR2(SF,i,j)-DR2(SF,i,j-1))*idy/hbar;         
         DR2(U,i,j) += pp;
      }
   }

   // get rid of SF memory now

   deallocate_field2D(SF);
}



void  init_initial_condition(initial_condition_info_t *IC, int nx, int ny)
{

   UU = allocate_field2D_grid(nx,ny, UGRID);
   VV = allocate_field2D_grid(nx+1,ny, VGRID);
   ETA = allocate_field2D_grid(nx-1,ny, ETAGRID);


   field2D_set_to_value(UU,0.0);
   field2D_set_to_value(VV,0.0);
   field2D_set_to_value(ETA,0.0);


  /* load in the u and v initial conditions */

  switch(IC->uic) {
  case 0:    break;                 // zero out things 
  case 1:    load_line_IC(IC->uic_file->str, UU);
             break;
  case 2:    load_full_IC(IC->uic_file->str, UU);
             break;
  default:   funwaveC_error_message("Bad u initial condition switch"); // put this in load_init_file
  }

   switch(IC->vic) {
   case -1:   break;                    // zero out things   
   case 0:    set_const_v(VV,atof(IC->vic_file->str)); 
              break;
   case 1:    load_line_IC(IC->vic_file->str, VV);
              break;
   case 2:    load_full_IC(IC->vic_file->str, VV);
              break;
   default:   funwaveC_error_message("Bad v initial condition switch");
   }

   switch(IC->etaic) {
   case -1:   break;                    // zero out things   
   case 0:    set_const_v(ETA, 0.0 ); 
              break;
   case 1:    load_line_IC(IC->etaic_file->str, ETA);
              break;
   case 2:    load_full_IC(IC->etaic_file->str, ETA);
              break;
   default:   funwaveC_error_message("Bad eta initial condition switch");
   }


   /* load the optional perturbation streamfunction */

   switch(IC->streamic) {
   case 0:    break;
   case 2:    load_streamfunction(IC->stream_file->str,UU,VV,H);
              break;
   default:   funwaveC_error_message("Wrong streamfunction flag in init file");
   }




}



void  deallocate_initial_condition()
{

  deallocate_field2D(UU);
  deallocate_field2D(VV);
  deallocate_field2D(ETA);


}