ezstep.h

这是一个生物的吸收仿真程序,在模拟方面也算的是一个不错的小程序.

/* ------------------------------------------------------------------------- 
 * ezstep.h -- Macros for EZSTEP
 *
 * Copyright (C) 1992, 1994, 1997, 1998 Dwight Barkley                 
 *
 * RCS Information
 * ---------------------------
 * $Revision: 3.1.1.1 $
 * $Date: 1998/06/02 21:07:53 $
 * ------------------------------------------------------------------------- */

#include "ezspiral.h"
#ifndef _EZSTEP_
#define _EZSTEP_



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

 * Define the macros U_PARTIAL_X(index) and U_PARTIAL_Y(index) for the 

 * computation of partial difference of U along X and Y axis respectively

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


#if U_PARTIAL_ON

#define U_PARTIAL_X(index)  dt * p1 * ((u[index] - u[index - I_INC]) + \
						(u[index + I_INC] - u[index])) / (2.0 * grid_h)

#define U_PARTIAL_Y(index)  dt * p2 * ((u[index] - u[index - J_INC]) + \
						(u[index + J_INC] - u[index])) / (2.0 * grid_h)

#endif



#if V_PARTIAL_ON

#define V_PARTIAL_X(index)  dt * p3 * ((v[index] - v[index - I_INC]) + \
						(v[index + I_INC] - v[index])) / (2.0 * grid_h)

#define V_PARTIAL_Y(index)  dt * p4 * ((v[index] - v[index - J_INC]) + \
						(v[index + J_INC] - v[index])) / (2.0 * grid_h)

#endif



/* ---------------------------------------------------------------------
 * Define the macros U_THRESHOLD(v) and G(u,v) according to the model
 * you wish to simulate. In principle these can be almost anything you
 * want.  Three examples are included.  
 * --------------------------------------------------------------------- */




#if 0  /* Standard Model */
  #define U_THRESHOLD(v)  ( one_o_a*(v) + b_o_a )
  #define G(u,v)          ( (u)-(v) )
#endif

#if 1  /* Simple Model for Turbulent Spirals */
  #define U_THRESHOLD(v)  ( one_o_a*(v) + b_o_a )
  #define G(u,v)          ( (u)*(u)*(u) - (v) )
#endif

#if 0  /* Bar and Eiswirth Model (Phys. Rev. E V. 48, p. 1635, 1993).
	* I do not include their case u>1 because it is not relevant. */
  #define U_THRESHOLD(v)  ( one_o_a*(v) + b_o_a )
  #define G(u,v)          (-(v) + ( ((u)<third) ? 0 : \
                                    (one-6.75*(u)*(1-(u))*(1-(u))) ) )  
#endif




/* ------------------------------------------------------------------------- 
 *          In general you should not need to change anything below          
 * ------------------------------------------------------------------------- */

/* I rely on the compiler to take care of most optimization.  The only place
 * where I give the compiler help is with the Laplacian formulas because
 * none of the compiles do the obvious thing.  */

/* ------------------------------------- 
 * Defining the kinetics macros:         
 * U_KINETICS(u,uth) and V_KINETICS(u,v) 
 * ------------------------------------- */

#if 1  
     /* Explicit u-kinetics */
  #define U_KINETICS(u, v ,uth) ( (u)+dt_o_eps*(u)*(one-(u))*((u)-(uth)) )

#endif

#if 0       
     /* Implicit u-kinetics */
     /* The original (Physica 49D) scheme can be obtained by defining
      * F2m and F2p to be one */
  #define F1m(u,uth) ( dt_o_eps*(one-(u))*((u)-(uth)) )
  #define F1p(u,uth) ( dt_o_eps*(u)*((u)-(uth)) )
  #define F2m(u,uth) ( one + dt_o_2eps*((uth)-(u)*(u)) )
  #define F2p(u,uth) ( one + dt_o_2eps*(two*(u) -(uth)-(u)*(u)) )

  #define U_KINETICS(u, v, uth) (                                      \
          ( (u) < (uth) ) ?                                        \
          (u)/(one-F1m(u,uth)*F2m(u,uth) ) :                       \
          ((u)+F1p(u,uth)*F2p(u,uth))/(one+F1p(u,uth)*F2p(u,uth)) )

#endif



#if 0

	#define PHI_T 0.01 + A * sin(omiga * dt * (Real)istep)

	#define U_KINETICS(u, v, uth)  ((u) + dt_o_eps * (u - u*u - (2*v + PHI_T) * (u - 0.002) / (u + 0.002)))

#endif

#define V_KINETICS(u,v) ( (v)+dt*G(u,v) )

/* ---------------------------------------- 
 * Defining the diffusion macros:           
 * U_DIFFUSION, V_DIFFUSION, ZERO_USED_SUMS 
 * ---------------------------------------- */

#if V_DIFF_ON        
     /* v is diffusing */
  #define U_DIFFUSION(index_1)     (dt_o_wh2 * Du * sigma_u[index_1]) 
  #define V_DIFFUSION(index_1)     (dtDv_o_wh2 * sigma_v[index_1]) 
  #define ZERO_USED_SUMS(index_1)  sigma_u[index_1] = sigma_v[index_1] = zero;
#else
     /* v is not diffusing */
  #define U_DIFFUSION(index_1)     (dt_o_wh2 * Du * sigma_u[index_1]) 
  #define V_DIFFUSION(index_1)     zero
  #define ZERO_USED_SUMS(index_1)  sigma_u[index_1] = zero;
#endif

/* -------------------------------- 
 * Defining the spatial-sum macros: 
 * ADD_TO_U_SUM, ADD_TO_V_SUM     
 * -------------------------------- */

#if NINEPOINT
/* 9-point Laplacian formula */
  #define ADD_TO_U_SUM(index,index_2) \
   {Real stupid_cc = u[index]; \
     sigma_u[(index_2)]  -= twenty*stupid_cc; \
     sigma_u[(index_2)+J_INC]  += four*stupid_cc; \
     sigma_u[(index_2)-J_INC]  += four*stupid_cc; \
     sigma_u[(index_2)+I_INC]  += four*stupid_cc; \
     sigma_u[(index_2)-I_INC]  += four*stupid_cc; \
     sigma_u[(index_2)+I_INC+J_INC] += stupid_cc; \
     sigma_u[(index_2)+I_INC-J_INC] += stupid_cc; \
     sigma_u[(index_2)-I_INC+J_INC] += stupid_cc; \
     sigma_u[(index_2)-I_INC-J_INC] += stupid_cc; \
   }
#else
/* 5-point Laplacian formula */
  #define ADD_TO_U_SUM(index,index_2) \
   {Real stupid_cc = u[index]; \
      sigma_u[(index_2)]  -= four*stupid_cc; \
      sigma_u[(index_2)+J_INC] += stupid_cc; \
      sigma_u[(index_2)-J_INC] += stupid_cc; \
      sigma_u[(index_2)+I_INC] += stupid_cc; \
      sigma_u[(index_2)-I_INC] += stupid_cc; \
   }
#endif         

#if !V_DIFF_ON   
/* v is not diffusing */
  #define ADD_TO_V_SUM(index,index_2) 
#else            
/* v is diffusing */
#if NINEPOINT
/* 9-point Laplacian formula */
  #define ADD_TO_V_SUM(index,index_2) \
   {Real stupid_cc = v[index]; \
     sigma_v[(index_2)]  -= twenty*stupid_cc; \
     sigma_v[(index_2)+J_INC]  += four*stupid_cc; \
     sigma_v[(index_2)-J_INC]  += four*stupid_cc; \
     sigma_v[(index_2)+I_INC]  += four*stupid_cc; \
     sigma_v[(index_2)-I_INC]  += four*stupid_cc; \
     sigma_v[(index_2)+I_INC+J_INC] += stupid_cc; \
     sigma_v[(index_2)+I_INC-J_INC] += stupid_cc; \
     sigma_v[(index_2)-I_INC+J_INC] += stupid_cc; \
     sigma_v[(index_2)-I_INC-J_INC] += stupid_cc; \
   }
#else          
/* 5-point Laplacian formula */
  #define ADD_TO_V_SUM(index,index_2) \
   {Real stupid_cc = v[index]; \
      sigma_v[(index_2)]  -= four*stupid_cc; \
      sigma_v[(index_2)+J_INC] += stupid_cc; \
      sigma_v[(index_2)-J_INC] += stupid_cc; \
      sigma_v[(index_2)+I_INC] += stupid_cc; \
      sigma_v[(index_2)-I_INC] += stupid_cc; \
  }
#endif 
#endif
#endif /*  _EZSTEP_  */