field.c

波浪数值模拟

      *dfd0++ = *fddn++;          /* then the near boundary */
    }
  }


}


void field2D_bary_v_to_eta(field2D *ABY, const field2D *A)
{
  int NY = ABY->N;
  int N = A->N;
  int M = A->M;
  const int MN = M-1;

  int i, j;

  const double *ad = &(A->data[0]);
  double *abyd = &(ABY->data[0]);

  double aminus, aplus;

  /* Check the grid types first - eventually put this into an #ifdef  */

  if ( (A->grid != VGRID) || (ABY->grid != ETAGRID) ) {
    field_error_message("Bad grid types to field2D_bary_v_to_eta()");
  }

  if (NY != N) {
    field_error_message("Bad N field2D sizes to field2D_bary_v_to_eta()");
  }

  if (ABY->M != M) {
    field_error_message("Bad M field2D sizes to field2D_bary_v_to_eta()");
  }

  for (i=0;i<N;i++) {
    aminus = *(ad+MN);
    aplus = *ad++;
    *abyd++ = 0.5 * (aminus+aplus);
    aminus = aplus; 
    for (j=1;j<M;j++) {
      aplus = *ad++;
      *abyd++ = 0.5 * (aplus + aminus);
      aminus = aplus;
    } 
  }


}

void field2D_bary_eta_to_v(field2D *ABY, const field2D *A)
{
  int NY = ABY->N;
  int N = A->N;
  int M = A->M;
  const int MN = M-1;
  int i, j;

  const double *ad = &(A->data[0]);
  double *abyd = &(ABY->data[0]);

  double aminus, aplus;

  /* Check the grid types first - eventually put this into an #ifdef  */

  if ( (A->grid != ETAGRID) || (ABY->grid != VGRID) ) {
    field_error_message("Bad grid types to field2D_bary_eta_to_v()");
  }

  if (NY != N) {
    field_error_message("Bad N field2D sizes to field2D_bary_eta_to_v()");
  }

  if (ABY->M != M) {
    field_error_message("Bad M field2D sizes to field2D_bary_eta_to_v()");
  }

  for (i=0;i<N;i++) {

    aminus = *ad++;
    for (j=0;j<M-1;j++) {
      aplus = *ad++;
      *abyd++ = 0.5 * (aplus + aminus);
      aminus = aplus;
    } 
    aplus = *(ad-M);
    *abyd++ = 0.5 * (aminus+aplus);

  }

}

/* This does the \bar^y() operation from a U point onto an VORT (vorticity) point */

void field2D_bary_u_to_vort(field2D *ABY, const field2D *A)
{
  int NY = ABY->N;
  int N = A->N;
  int M = A->M;
  const int MN = M-1;
  int i, j;

  const double *ad = &(A->data[0]);
  double *abyd = &(ABY->data[0]);

  double aminus, aplus;

  /* Check the grid types first - eventually put this into an #ifdef  */

  if ( (A->grid != UGRID) || (ABY->grid != VORTGRID) ) {
    field_error_message("Bad grid types to field2D_bary_u_to_vort()");
  }

  if (NY != N) {
    field_error_message("Bad N field2D sizes to field2D_bary_u_to_vort()");
  }

  if (ABY->M != M) {
    field_error_message("Bad M field2D sizes to field2D_bary_u_to_vort()");
  }

  for (i=0;i<N;i++) {

    aminus = *ad++;
    for (j=0;j<M-1;j++) {
      aplus = *ad++;
      *abyd++ = 0.5 * (aplus + aminus);
      aminus = aplus;
    } 
    aplus = *(ad-M);
    *abyd++ = 0.5 * (aminus+aplus);

  }

}

/* Puts a ETA-grid  field2D onto a VORT-grid field2D 
    Calculates a \bar^{xy} A and puts it into ABXY 
   - requires ABXY have size (N-1,M) or relative to A (N,M)   */

void field2D_barxy_eta_to_vort(field2D *ABXY, const field2D *A)
{
  int NXY = ABXY->N;
  int N = A->N;
  int M = A->M;
  //  const int MN = M-1;
  int i, j;

  const double *adup = &(A->data[M]);
  const double *addn = &(A->data[0]);

  double *abxyd = &(ABXY->data[0]);

  double abxp, abxm; 

  /* Check the grid types first - eventually put this into an #ifdef  */

  if ( (A->grid != ETAGRID) || (ABXY->grid != VORTGRID) ) {
    field_error_message("Bad grid types to field2D_barxy_eta_to_vort()");
  }

  if (ABXY->M != M) {
    field_error_message("Bad M field2D sizes to field2D_barxy_eta_to_vort()");
  }

  switch (NXY - N)  {
  case 1: abxyd = REF2(ABXY,1);
    break;
  case -1: abxyd = REF2(ABXY,0);
    break;
  default:   field_error_message("Bad N field2D sizes to field2D_barxy_eta_to_vort()");
  }

  /* Recall that vort points are on the same j grid as V points so the wrap around comes at the end */

  for (i=0;i<N-1;i++) {
    abxm = ( *adup++ + *addn++);          /* first calculate \bar^x(A) at j=0 */
    for (j=0;j<M-1;j++) { 
      abxp = ( *adup++ + *addn++);        /* calculate \bar^x(A) at j+1 */
      *abxyd++ = 0.25 * (abxm + abxp);    /* average \bar^y ( \bar^x(A)) at j & j+1 */
      abxm = abxp;                        /* save j+1 value into j value */
    }
    abxp = ( *(adup-M) + *(addn-M));      /* take care of the boundary */
    *abxyd++ = 0.25 * (abxm + abxp);
  }



}



/* Calculates a \delta_x difference of F and puts it into DF 
   - requires DF have size (N-1,M) or (N+1,M) relative to F (N,M)        
   - if size=(N+1,M) then the i=0 & N points are zeros!           */

void field2D_diffx(field2D *DF, const field2D *F, double idx)
{
  int NF = F->N;
  int NDF = DF->N;
  const int M = F->M;

  int i, j;

  const double *fdup = &(F->data[M]);
  const double *fddn = &(F->data[0]);
  double *dfd = &(DF->data[0]);
  double *dfd0, *dfdN;

  if ( (NDF != NF-1) && (NDF != NF+1) ) {
    field_error_message("Bad nx field2D sizes to field2D_diffx()");
  }

  if (DF->M != M) {
    field_error_message("Bad M field2D sizes to field2D_diffx()");
  }

  if (NDF == NF-1) {
    dfd = &(DF->data[0]);
    for (i=0;i<NF-1;i++) {
      for (j=0;j<M;j++) {
        *dfd++ = idx * ( *fdup++ - *fddn++);
      }
    }
  }
  else {
    dfd = &(DF->data[M]);
    for (i=0;i<NF-1;i++) {
      for (j=0;j<M;j++) {
        *dfd++ = idx * ( *fdup++ - *fddn++);
      }
    }
    dfd0 = DF->row[0];
    dfdN = DF->row[NDF-1];
    for (j=0; j<M; j++) {
      *dfd0++ = 0.0;
      *dfdN++ = 0.0;
    }
  }

}



/* Calculates a \delta_x difference of F and puts it into DF 
   - but uses a 4th order scheme   [1 -13.5 13.5  -1]/(21/2) away from the boundaries
   - the error here should be O(dx^4)
   - near the boundaries it's still 2nd order.
   - requires DF have size (N-1,M) or (N+1,M) relative to F (N,M)        
   - if size=(N+1,M) then the i=0 & N points are zeros!           */

void field2D_diffx_4th_order(field2D *DF, const field2D *F, double idx)
{
  int NF = F->N;
  int NDF = DF->N;
  const int M = F->M;

  int i, j;

  const double *fdup = &(F->data[M]);
  const double *fddn = &(F->data[0]);
  double *dfd = &(DF->data[0]);

  const double ifac = 2.0/21.0;

  if ( (NDF != NF-1) && (NDF != NF+1) ) {
    field_error_message("Bad nx field2D sizes to field2D_diffx()");
  }

  if (DF->M != M) {
    field_error_message("Bad M field2D sizes to field2D_diffx()");
  }

  if (NDF == NF-1) {
    dfd = &(DF->data[0]);
  }
  else {
    dfd = &(DF->data[M]);
  }


  fprintf(stderr,"4th order dx differencing not yet tested!");

  /* First do the near x boundary points in same 2nd order fashion */


  i = 0;
  for (j=0;j<M;j++) {
      *dfd++ = idx * ( *fdup++ - *fddn++);
    }

  for (i=1;i<NF-2;i++) {
    for (j=0;j<M;j++) {
      *dfd++ = idx * ifac * ( 13.5 * ( *fdup - *fddn) - (*(fdup+M) - *(fddn-M) ));
      fdup++;
      fddn++;
    }
  }

  i = NF-2;   /* this si the far boundary - 2nd order here */
  for (j=0;j<M;j++) {
      *dfd++ = idx * ( *fdup++ - *fddn++);
    }


}



/* Calculates  DF =  \delta_x \delta_x (F) 
   - requires DF have size (N-2,M) relative to F (N,M)        */


void field2D_diffx2(field2D *FXX, const field2D *F, double idx)
{
  int NF = F->N;
  int NFXX = FXX->N;
  const int M = F->M;

  int i, j;

  const double *fdup = &(F->data[2*M]);
  const double *fdmid = &(F->data[M]);
  const double *fddn = &(F->data[0]);

  double *fxxd = &(FXX->data[0]);

  const double idx2 = idx*idx;

  if ( (NFXX != NF-2)  ) {
    field_error_message("Bad nx field2D sizes to field2D_diffx2()");
  }

  if (FXX->M != M) {
    field_error_message("Bad M field2D sizes to field2D_diffx2()");
  }

  for (i=0;i<NF-2;i++) {
    for (j=0;j<M;j++) {
      *fxxd++ = idx2 * ( *fdup++ - 2.0*( *fdmid++) + *fddn++);
    }
  }

}

/* This calculates  DF += \delta_x(F)  but requires that the DF(N-1,M) when F(N,M)  */

void field2D_diffx_add(field2D *DF, const field2D *F, double idx)
{
  int NXf = F->N;
  int NXdf = DF->N;
  const int M = F->M;

  int i, j;

  const double *fdup = &(F->data[M]);
  const double *fddn = &(F->data[0]);
  double *dfd = &(DF->data[0]);

  if (NXdf != NXf-1) {
    field_error_message("Bad nx field2D sizes to field2D_diffx()");
  }

  if (DF->M != M) {
    field_error_message("Bad M field2D sizes to field2D_diffx()");
  }

  for (i=0;i<NXf-1;i++) {
    for (j=0;j<M;j++) {
      *dfd++ += idx * ( *fdup++ - *fddn++);
    }

  }
}


/* Calculates a \delta_x difference of F and puts it into DF, converts from eta points to V points
   - requires DF have size (N,M) relative to F (N,M)        */

void field2D_diffy_eta_to_v(field2D *DF, const field2D *F, double idy)
{
  int N = F->N;
  int Nf = DF->N;
  int M = F->M;
  const int MN = M-1;

  int i, j;

  const double *fd = &(F->data[0]);
  double *dfd = &(DF->data[0]);
  double fplus, fminus;

  /* Check the grid types first - eventually put this into an #ifdef  */

  if ( !((DF->grid == VGRID)||(DF->grid == VTGRID)) || (F->grid != ETAGRID) ) {
    field_error_message("Bad grid types to field2D_diffy_eta_to_v()");
  }

  if (N != Nf) {
    field_error_message("Bad nx field2D sizes to field2D_diffy_eta_to_v()");
  }

  if (DF->M != M) {
    field_error_message("Bad M field2D sizes to field2D_diffy_eta_to_v()");
  }

  for (i=0;i<N;i++) {
  
    fminus = *fd++;
    for (j=0;j<M-1;j++) {
      fplus = *fd++;
      *dfd++ = idy * ( fplus - fminus );
      fminus = fplus;
    }
    fplus = *(fd-M);
    *dfd++ = idy * ( fplus - fminus );
  }

  /*  for (i=0;i<N;i++) {

    for (j=0;j<M-1;j++) {
      *dfd++ = idy * ( *(fd+1) - *fd);
      fd++;
    }
    *dfd++ = idy * ( *(fd-MN) - *fd);
    fd++;
    }  */
}


void field2D_diffy_eta_to_v_add(field2D *DF, const field2D *F, double idy)
{
  int N = F->N;
  int Nf = DF->N;
  int M = F->M;
  const int MN = M-1;

  int i, j;

  const double *fd = &(F->data[0]);
  double *dfd = &(DF->data[0]);
  double fminus, fplus;

  /* Check the grid types first - eventually put this into an #ifdef  */

  if ( (F->grid != ETAGRID) || (DF->grid != VGRID) ) {
    field_error_message("Bad grid types to field2D_diffy_eta_to_v_add()");
  }

  if (N != Nf) {
    field_error_message("Bad nx field2D sizes to field2D_diffy_eta_to_v_add()");
  }

  if (DF->M != M) {
    field_error_message("Bad M field2D sizes to field2D_diffy_eta_to_v_add()");
  }

  for (i=0;i<N;i++) {
  
    fminus = *fd++;
    for (j=0;j<M-1;j++) {
      fplus = *fd++;
      *dfd++ += idy * ( fplus - fminus );
      fminus = fplus;
    }
    fplus = *(fd-M);
    *dfd++ += idy * ( fplus - fminus );
  }

  /*
  for (i=0;i<N;i++) {

    for (j=0;j<M-1;j++) {
      *dfd++ += idy * ( *(fd+1) - *fd);
      fd++;
    }
    *dfd++ += idy * ( *(fd-MN) - *fd);
    fd++;
  }
  */
}


/* Calculates a \delta_x difference of F and puts it into DF, converts from eta points to V points
   - requires DF have size (N,M) relative to F (N,M)        */

void field2D_diffy_v_to_eta(field2D *DF, const field2D *F, const double idy)
{
  const int NF = F->N;
  const int NDF = DF->N;
  const int M = F->M;
  const int MN = M-1;

  int i, j;

  const double *fd;
  double *dfd = &(DF->data[0]);
  double fplus, fminus;

  /* Check the grid types first - eventually put this into an #ifdef  */

  if ( ( (F->grid != VGRID)&&(F->grid !=VTGRID) ) || (DF->grid != ETAGRID) ) {
    field_error_message("Bad grid types to field2D_diffy_v_to_eta()");
  }

  if (DF->M != M) {
    field_error_message("Bad M field2D sizes to field2D_diffy_v_to_eta()");
  }

  switch (NF-NDF) {
  case 0: fd = REF2(F, 0);
    break;
  case 2: fd = REF2(F, 1);
    break;
  default:  field_error_message("Bad nx field2D sizes to field2D_diffy_v_to_eta()");
  }


  for (i=0;i<NDF;i++) {
    fminus = *(fd+MN);
    fplus = *fd++;
    *dfd++ = idy * ( fplus - fminus );
    fminus = fplus;

    for (j=1;j<M;j++) {
      fplus = *fd++;
      *dfd++ = idy * ( fplus - fminus );
      fminus = fplus;
    }
  }
}


void field2D_diffy_v_to_eta_add(field2D *DF, const field2D *F, const double idy)
{
  const int NF = F->N;
  const int NDF = DF->N;
  const int M = F->M;
  const int MN = M-1;

  int i, j;

  const double *fd;
  double *dfd = &(DF->data[0]);
  double fplus, fminus;

  /* Check the grid types first - eventually put this into an #ifdef  */

  if ( !((F->grid == VGRID)||(F->grid == VTGRID)) || (DF->grid != ETAGRID) ) {
    field_error_message("Bad grid types to field2D_diffy_v_to_eta_add()");
  }

  if (DF->M != M) {
    field_error_message("Bad M field2D sizes to field2D_diffy_eta_to_v_add()");
  }

  switch (NF-NDF) {
  case 0: fd = REF2(F, 0);
    break;
  case 2: fd = REF2(F, 1);
    break;
  default:  field_error_message("Bad nx field2D sizes to field2D_diffy_v_to_eta()");
  }

  for (i=0;i<NDF;i++) {
    fminus = *(fd+MN);
    fplus = *fd++;
    *dfd++ += idy * ( fplus - fminus );
    fminus = fplus;

    for (j=1;j<M;j++) {
      fplus = *fd++;
      *dfd++ += idy * ( fplus - fminus );
      fminus = fplus;
    }
  }