routines.c

来自「波浪数值模拟」· C语言 代码 · 共 161 行

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

/* ---- routines.c   */

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

#include "routines.h"
//#include "misc.h"
//#include "dynamics.h"



double total_kinetic_energy()
{
  int i,j;
  double tot=0.0;

  /* Here I should really multiply by the water depth to get the
     depth integrated kinetic energy.......        */

  for (i=1;i<nx;i++)
    for (j=0;j<ny;j++)
      tot += kinetic_energy(i,j,UU,VV,H);

  tot *= dx*dy;

  return tot;

}


double total_potential_energy()
{
  int i,j;
  double tot=0.0;

  /* Here I should really multiply by the water depth to get the
     depth integrated kinetic energy.......        */

  for (i=1;i<nx;i++)
    for (j=0;j<ny;j++)
      tot += potential_energy(i,j,ETA);

  tot *= dx*dy;

  return tot;

}


double kinetic_energy(const int i, const int j, const field2D *U, const field2D *V, const field2D *H)
{
    double ubar,vbar,tot;
    int b0,b1;

    switch (j) {
      case 0:    b0=ny-1;
                 b1=0;
                 break;
      default:   b0=j-1;
                 b1=j;
      }
    
    vbar = 0.5*(DR2(V,i,b0)+DR2(V,i,b1));
    ubar = 0.5*(DR2(U,i-1,j)+DR2(U,i,j));

    tot = 0.5* DR2(H,i,j) * (SQR(vbar)+SQR(ubar));

    return tot;

}


double potential_energy(int i,int j, const field2D *ETA)
{
    double tot;
    double eta;
    double g = 9.81;

    eta = DR2(ETA,i-1,j);
    tot = 0.5*g*eta*eta;

    return tot;
}    






field2D*  calculate_vorticity(const field2D *U, const field2D *V)
{

  int i, j, nx, ny;
  static field2D *VORT = NULL;
  double *vortd;
  int jp;

  nx = V->N -1;
  ny = V->M; 

  /* potential memory leak here until program ends!! */

  if ( NOT_ALLOCATED( VORT ) ) {
    VORT = allocate_field2D(nx,ny);  
  }

  vortd = VORT->data;

  for (i=0;i<nx;i++) {
    for (j=0;j<ny;j++) {
       jp = JPM(j);
       *vortd++ = idx*(DR2(V,i+1,j)-DR2(V,i,j)) - idy*(DR2(U,i,jp)-DR2(U,i,j));
    }
  }

  return VORT;
        
}


double vorticity(const int i, const int j, const field2D *U, const field2D *V)
{

    double t0;
    int jp;

    jp = JPM(j);

    /*
    switch(j) {
    case  ny-1:  jp=0;  break;
    default: jp=j+1; 
    }
    */
    t0 = idx*(DR2(V,i+1,j)-DR2(V,i,j)) - idy*(DR2(U,i,jp)-DR2(U,i,j));
    return t0;
        
}