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


/*    
      funwaveC.c  -  nonlinear bousinessq wave model based on Jim Kirby et al.'s funwave

*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <pthread.h>
#include <assert.h>

#include "bdefs.h"
#include "cmdline.h"
#include "gui.h"
#include "input.h"
#include "flags.h"
#include "funwaveC_timestep.h"
#include "depth.h"
#include "floats.h"
#include "depth.h"
#include "eta_source.h"
#include "sponge.h"
#include "lateral_mixing.h"
#include "bottom_stress.h"
#include "forcing.h"             /* this is the external forcing such as wind or rad stress grad if appropriate */
#include "breaking.h"
#include "initial_condition.h"
#include "routines.h"
#include "output.h"


/* Here the global model domain size parameters are defined */

int nx;        /* number of U x-grid points, ETA is nx-1, V is nx+1 */
int ny;        /* number of (U,V,ETA,...) y-grid points */

double dx;     /* x-grid spacing units of m */
double dy;     /* y-grid spacing units of m */

double idx;    /* =1/dx: This is in units of 1/m, grid spacing  */
double idy;    /* =1/dy:                                        */
double idx2;   /* = idx*idx  - used for 2nd derivatives */
double idy2;   /* = idy*idy  - used for 2nd derivatives */



const double f = 1.0e-4;             // Coriolis 
const double gravity = 9.81;         // m/s^2 

double global_time = 0.0;

double nu_newt;
double nu_bi;



field2D *UU;
field2D *VV;
field2D *ETA;

field2D *FXX, *FYY;      // x and y  depth-integrated forcing


/*  Variables for bottom stress */

field2D *TAUBX, *TAUBY;

field2D *H;      /* This is the bathymetry */

/* Explicit declaration of funwaveC exit - need to exit this way to make sure GUI window is cleaned up */

void funwaveC_exit(int exitcode);

void funwaveC_error_message(char *s)
{

  fprintf(stderr,"** funwaveC: ERROR: %s\n",s);
  funwaveC_exit(-1);

}

void funwaveC_error_message_file_line(char *s, char *file, int lineno)
{

  fprintf(stderr,"** funwaveC: ERROR: in file %s, line %d \n: %s\n",file,lineno,s);
  funwaveC_exit(-1);

}

void funwaveC_exit_message(char *s)
{

  fprintf(stderr,"funwaveC exit:  %s\n",s);
  funwaveC_exit(0);

}




void  deallocate_all_fields()
{

  deallocate_initial_condition();
  deallocate_depth_variables();
  deallocate_funwaveC_timestep();        // free the memory with these arrays
  deallocate_lateral_mixing();

  deallocate_passive_tracer_variables();
  free_floats();



}


void initialize_all(inputinfo *I)
{

  /* first set up the global variables includes 
     grid sizes (nx,ny), grid space (dx, dy), and inverse grid space (idx, idy) */

  nx = I->nx;   ny = I->ny;

  /* check to make sure that ny is large enough  */

  if ((ny < 6)) {
       funwaveC_error_message("* ny is too small !\n");
       }

  dx = I->dx;   dy = I->dy;
  idx = 1.0/dx;   idy = 1.0/dy;
  idx2 = idx*idx; idy2 = idy*idy;


  /* Now initialize the bathymetry file: this routine is in depth.c  */

  init_depth_variables(I->D, nx, ny);

  /* Next initialize the U, V, and ETA fields */

  init_initial_condition(I->IC, nx, ny);

  /* Next initialize all the funwaveC time stepping stuff */

  init_funwaveC_timestep(nx, ny, I->Dynamics);

  /* The wave maker: eta_source */

  init_eta_source(I->ES, nx, ny, I->T->dt);

  /* Initialize the sponge layer */

  init_sponge_layer(I->SP, nx, ny);

  /* Initialize the wave breaking terms */

  init_breaking(I->BR, Depth_Vars);

  /* next initialize the bottom stress, lateral_mixing and forcing */

  c_d = I->cd;
  init_bottom_stress( I->cd) ;

  /* Initialize the lateral mixing terms - newtonian or biharmonic */

  init_lateral_mixing(nx, ny, I->DT, I->nu_newt, I->nu_bi);

  /* Now initialize the forcing - either a wind stress or other body force */
  /* needs to have bottom stress, depth, and IC initialized */

  init_forcing(I->F, nx, ny);  

  /* lastly initialize the passive tracer and floats  */

  init_passive_tracer_variables(nx,ny, I->TR);

  init_floats(I->FL);

}


/* ---------------------------------------------------------------- */


// REPORTS on the definitions turned on

void  defs_report(dynamics_t Dynamics)
{

    // First report on time-stepping and dx==dy diagnostics

    fprintf(stderr,"-------------------------------\n");
    fprintf(stderr,"* Time-stepping with AB3.\n");

    switch (Dynamics) {
    case LINEAR_DYNAMICS: 
      fprintf(stderr,"* Linear Shallow Water Eq Dynamics \n");
      break;

    case NSWE_DYNAMICS:
      fprintf(stderr,"* Nonlinear Shallow Water Eq Dynamics \n");
      break;

    case PEREGRINE_DYNAMICS:
      fprintf(stderr,"* Linear Peregrine Bousinessq Dynamics \n");
      break;

    case NWOGU_DYNAMICS:
      fprintf(stderr,"* Linear Nwogu Bousinessq  Dynamics \n");
      break;

    case WEI_KIRBY_DYNAMICS:
      fprintf(stderr,"* Fully nonlinear Wei & Kirby Bousinessq  Dynamics \n");
      break;
    }


    fprintf(stderr,"-------------------------------\n");

    
}


/* (1) Routine to report on the diagnostics of the model run  */

void stability_report(double dt, double dx, double dy, double cd, field2D *U, field2D *V, field2D *H)
{

    double tmp, t1, Lx, Ly, L, Umax;
    double CFL_UVx, CFL_WAVEx;
    double CFL_UVy, CFL_WAVEy;
    double uu,vv;
    double min_depth, min_h0;
    int i,j;

    // compute the Lx and  Ly domain sizes

    Lx = (nx-1.0)*dx;
    Ly = ny*dy;

    if (Lx>Ly) L = Ly; else L=Lx;

    // find the max U current and smallest depth, do this ad hoc
    
    tmp=0.0;
    for (i=1;i<nx-1;i++)
      for (j=0;j<ny;j++) {
	uu = DR2(U,i,j);   vv = DR2(V,i,j);
          t1  = uu*uu + vv*vv;
          if (t1>tmp)  tmp=t1;
      }
    Umax = sqrt(tmp);

    min_h0 = min_field2D(H);

    /* Do the CFL #  - this should be < 0.24 (AB2)  < 0.2 (AB3)  */

    CFL_UVx = Umax*dt*idx;
    CFL_WAVEx = sqrt(gravity*max_field2D(H)) * dt*idx;
    CFL_UVy = Umax*dt*idy;
    CFL_WAVEy = sqrt(gravity*max_field2D(H)) * dt*idy;

    fprintf(stderr,"-----------* Stability Report *-------------------------------------------\n");
    fprintf(stderr,"Max U = %g (m/s),  dt=%g (sec),  dx=%g (m) Lx=%g (m) Ly=%g (m)\n",Umax,dt,dx,Lx,Ly);
    fprintf(stderr,"cd=%g,   nu_bi=%g (m^4/s), nu_newt=%g (m^2/s), min h=%g (m)\n",cd,nu_bi,nu_newt,min_h0);


    /* Here the models Courant-Freidrichs-Levy numbers for both advection & shallow water
       waves are reported on (for x and y) and warnings given                  */

    fprintf(stderr,"CFL_UV-x # = %g , CFL_UV-y # = %g :",CFL_UVx, CFL_UVy);

    if ( (CFL_UVx > 0.2) || (CFL_UVy > 0.2) )
       fprintf(stderr,"Warning: model may be CFL UV unstable!\n");
    else
       fprintf(stderr,"CFL UV number OK\n");

    fprintf(stderr,"CFL_WAVEx # = %g , CFL_WAVEy # = %g  :",CFL_WAVEx, CFL_WAVEy);

    if ( (CFL_WAVEx > 0.4) || (CFL_WAVEy > 0.4) ) {
      funwaveC_error_message("** CFL WAVE NUMBER is too large. QUIT");
    }

    if ( (CFL_WAVEx > 0.2) || (CFL_WAVEy > 0.2) )
       fprintf(stderr,"Warning: model may be CFL WAVE unstable!\n");
    else
       fprintf(stderr,"CFL WAVE number OK\n");

    /* Now check the water depths */

    min_depth = check_depth(ETA);
    fprintf(stderr,"* Min depth (h+eta) = %g (m)\n", min_depth);
    report_mass_conservation(ETA);

    lateral_mixing_stability_report( Umax, L, dt );   /* report on lateral mixing stability parameters */

    bottom_stress_stability_report( Umax, min_h0, dt );   /* report on lateral mixing stability parameters */

    sponge_layer_stability_report(dt);  /* Report on Sponge Layer Stability */

    breaking_report();    /* Report on Wave breaking - whether it happenned or not */

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


//void diagnostic_report( )


/* --------- Integration Routines ------------------------------------ */




void funwaveC_level1_integrate(timing *T, cmdline *Cmd, oflags *level1)
{
  int i;

  for (i=0;i<T->level1_iter;i++) {
     funwaveC_timestep_AB3(T);
    }                       //  end of level 1 loop:  for(i=0.....  

  global_time = T->current_time;

  if (Cmd->verbose){
    verbose_output(T, total_kinetic_energy() , total_potential_energy() );
  }

  level_output(level1);      // level1 output stuff

  if (Cmd->diagnostic==1) {
    field2D_memory_usage();
    stability_report(T->dt,dx,dy,c_d,UU,VV,H);   
  }

}


void funwaveC_level2_integrate(timing *T, cmdline *Cmd, outputflags *A)
{
  int k;

  for (k=0;k<T->level2_iter;k++) {

       funwaveC_level1_integrate(T,Cmd,A->F1);

     }                               //  end of level 2 loop: for(k=0......

  level_output(A->F2);      // level2 output stuff

  if (Cmd->diagnostic==2) {
    field2D_memory_usage();
    stability_report(T->dt,dx,dy,c_d,UU,VV,H);   // final stability report
  }

}


void funwaveC_level3_integrate(timing *T, cmdline *Cmd, outputflags *A)
{
  int j;

  funwaveC_timestep_euler(T);      // do the first time step 

  global_time = T->current_time;
  level_output(A->F1);      // level1 output stuff

  if (Cmd->verbose) 
    verbose_output(T,total_kinetic_energy(),  total_potential_energy());

  funwaveC_timestep_AB2(T);  // do the next two time steps
  global_time = T->current_time;
  //  level_output(A->F1);      // level1 output stuff

  if (Cmd->verbose) 
    verbose_output(T,total_kinetic_energy(),  total_potential_energy());


  /* iterate through the rest of the time steps write kinetic energy if verbose*/
  /* This is the heart of the code */

  for (j=0;j<T->level3_iter;j++) {

    funwaveC_level2_integrate(T, Cmd, A);

    }                                  // end of loop: for(j=0....

  level_output(A->F3);

  if (Cmd->diagnostic==3) {
    field2D_memory_usage();
    stability_report(T->dt,dx,dy,c_d,UU,VV,H);   // final stability report
  }

}



#if THREAD_YES==1
pthread_t gui_thread;
pthread_t main_thread;
void *gui_thread_result;
switch_t gui_flag = OFF;
#endif /* THREAD_YES */


/* Generic exit function that cleans up the various threads before 
   exiting main.    This doesn't work yet           */

void funwaveC_exit(int exitcode)
{

#if THREAD_YES==1

  if (gui_flag == ON) {
    pthread_cancel(gui_thread);     // first the thread must be signalled to die 
    pthread_join(gui_thread,&gui_thread_result);     // Then we have to collect the thread 
  }

#endif /* THREAD_YES */

  exit(exitcode);

}



int main(int argc,char *argv[])
{
  gint status;

  inputinfo *I;
  outputflags *A;
  timing *T;
  cmdline *Cmd;


  fprintf(stderr,"******* funwaveC  %d.%d.%d *******************\n",
          FUNWAVEC_MAJOR_VERSION_NUMBER, FUNWAVEC_MINOR_VERSION_NUMBER, FUNWAVEC_MICRO_VERSION_NUMBER);

#if THREAD_YES
  gtk_init(&argc,&argv);

  main_thread = pthread_self();
#endif /*THREAD YES */

  Cmd = process_args2(argc,argv);

  I = load_init_file(Cmd->init_file);
  A = I->A;
  T = I->T;

  if (Cmd->parse) {
    funwaveC_exit_message("Finished Parsing init file\n");
  }

#if THREAD_YES==0
  if (Cmd->gui) {
    fprintf(stderr,"** Warning: No GUI output possible w/ no thread capability compiled in!\n");
  }
#endif /* THREAD_YES */


  /* Initialize Everything: depth, IC, dynamics, floats, tracer etc. */

  initialize_all(I);


  if (Cmd->diagnostic) {
    defs_report( I->Dynamics );           // report on what model options are set
    field2D_memory_usage();              // roughly estimate memory requirements  
    timing_report(T);              // report on the timing variables
    stability_report(T->dt,dx,dy,c_d,UU,VV,H);   // check the stability of the initial conditions
    eta_source_report();             /* Report on the wave maker stuff */
    passive_tracer_report(I->TR);  // report on the passive tracer variables
    floats_report(I->FL);          // report on the floats in the model 
  }


  /* Start the Gui timer */
#if THREAD_YES==1
  if (Cmd->gui) {
    T->msg = Cmd->msg;
    T->message = Cmd->message;
    T->initfilename = Cmd->init_file;
    gui_flag = ON;
    status = pthread_create(&gui_thread, NULL, (void *) gui_start, T);
  }
#endif /* THREAD_YES */


  /* do the integration */

  funwaveC_level3_integrate(T,Cmd,A);


  /* free all the allocated memory */

 deallocate_all_fields();   // deallocate all the field2D stuff
 outputflags_destroy(A);    // free the oflags memory


 funwaveC_exit(0);   // We exit here by cleaning up the GUI thread if necessary

 return 0;

}