C++ vectorField::lbound方法代码示例

/*
 * Adjust the time step according to the CFL stability criterion
 */
void adjustTimeStep(vectorField& V)
{
    // Find maximum velocity magnitude
    double maxV = 0.0;

    // NEEDS_WORK: Blitz should provide a norm(vectorField) function.
    // This is ugly.

    for (int i=V.lbound(0); i <= V.ubound(0); ++i)
      for (int j=V.lbound(1); j <= V.ubound(1); ++j)
        for (int k=V.lbound(2); k <= V.ubound(2); ++k)
        {
            double normV = norm(V(i,j,k));
            if (normV > maxV)
                maxV = normV;
        }

    cout << "Maximum velocity is " << maxV << " m/s" << endl;

    maxV += 1e-10;   // Avoid divide-by-zero

    // Steve K: need to have spatialStep^2 / diffusion constant
    // diffusion constant = eta * recip_rho

    delta_t = 0.3 * spatialStep / maxV;
    const double maxTimeStep = 0.01;

    if (delta_t > maxTimeStep)
        delta_t = maxTimeStep;

    cout << "Set time step to " << delta_t << " s" << endl;
}

BZ_END_STENCIL


/*
 * Allocate arrays and set their initial state
 */
void setup(const int N, vectorField& V, vectorField& nextV, scalarField& P,
    scalarField& P_rhs, vectorField& advect, vectorField& force)
{
    // A 1m x 1m x 1m domain
    spatialStep = 1.0 / (N - 1);
    geom = UniformCubicGeometry<3>(spatialStep);

    // Allocate arrays
    allocateArrays(shape(N,N,N), advect, V, nextV, force);  // vector fields
    allocateArrays(shape(N,N,N), P, P_rhs);                 // scalar fields

    // Since incompressibility is assumed, pressure only shows up as
    // derivative terms in the equations.  We choose airPressure = 0
    // as an arbitrary datum.

    airPressure = 0;             // Pa
    rho = 1000;                  // density of fluid, kg/m^3
    recip_rho = 1.0 / rho;       // inverse of density
    eta = 1.0e-6;                // kinematic viscosity of fluid, m^2/s
    gravity = 9.81;              // m/s^2
    delta_t = 0.001;             // initial time step, in seconds
    volume = pow3(spatialStep);  // cubic volume associated with grid point

    // Kludge: Set eta high, so that the flow will spread faster.
    // This means the cube is filled with molasses, rather than water.
    eta *= 1000;
   
    // Initial conditions: quiescent
    V = 0.0; 
    P_rhs = 0.0;
    advect = 0.0;
    nextV = 0.0;

    // Initial pressure condition: gravity causes a linear increase
    // in pressure with depth.  Note that tensor::k means the index
    // associated with the z axis (they are labelled i, j, k).
#ifdef NO_GRAVITY
    gravityPressureGradient = 0.0;
    P = 0.0;
#else
    gravityPressureGradient = spatialStep * gravity * rho;

#ifdef BZ_HAVE_NAMESPACES
    P = airPressure + tensor::k * gravityPressureGradient;
#else
    P = airPressure + k * gravityPressureGradient;
#endif

#endif

    // Set up the forcing function: gravity plus a stirring force
    // at the bottom
    double gravity_z = gravity * rho;

    const int x = 0, y = 1, z = 2;
    force[x] = 0.0;
    force[y] = 0.0;
#ifdef NO_GRAVITY
    force[z] = 0.0;
#else
    force[z] = gravity_z;    
#endif

#ifndef NO_STIRRING
    // Centre of the stirring
    int centrex = int(2 * N / 3.0);
    int centrey = int(2 * N / 3.0);
    int centrez = int(4 * N / 5.0);

    const double stirRadius = 1.0 / 3.0;

    vector3d zaxis(0,0,1);

    // Loop through the 2D slice where the stirring occurs

    for (int i=force.lbound(firstDim); i <= force.ubound(firstDim); ++i)
    {
      for (int j=force.lbound(secondDim); j <= force.ubound(secondDim); ++j)
      {
         // Vector from the centre of the stirring to the current
         // coordinate

         vector3d r((i-centrex) * spatialStep, (j-centrey) * spatialStep, 0.0);

         if (norm(r) < stirRadius)
         {
             // The cross product of the z-axis and the vector3d to this
             // coordinate yields the direction of the force.  Multiply
             // by gravity to get a reasonable magnitude (max force =
             // 5 * gravity)
             force(i,j,centrez) += cross(zaxis, r) 
                 * (5 * gravity_z / stirRadius);
         }
      }
//.........这里部分代码省略.........