field.c

波浪数值模拟

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

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




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

void field2D_diffy_u_to_vort(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 ( (F->grid != UGRID) || (DF->grid != VORTGRID) ) {
    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_u_to_vort()");
  }

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

  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 VORT points to U points
   - requires DF have size (N,M) or (N-2,M) relative to F (N,M)        */

void field2D_diffy_vort_to_u(field2D *DF, const field2D *F, const double idy)
{
  const int NF = F->N;
  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 fminus, fplus;

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

  if ( (F->grid != VORTGRID) || !( (DF->grid == UGRID)||(DF->grid == UTGRID)) ) {
    field_error_message("Bad grid types to field2D_diffy_vort_to_u()");
  }

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

  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_vort_to_u()");
  }

  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;
      }
  }
  /*
  for (i=0;i<NDF;i++) {
    *dfd++ = idy * ( *fd - *(fd+MN));
    fd++;

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



/* This function takes input F and calculates  \delta_x(F) and \delta_y(F) in one go to make things a tad faster */

void field2D_gradient_eta_to_uv(field2D *FX, field2D *FY, const field2D *F, const double idx, const double idy)
{
  int N = F->N;
  int NFY = FY->N;
  int NFX = FX->N;
  int M = F->M;

  const int MN = M-1;
  double fplus, fminus, fup;
  int i, j;

  const double *fd = &(F->data[0]);
  double *fyd = REF2(FY,0);
  double *fxd = REF2(FX,0);

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

  if ( !((FX->grid == UTGRID)||(FY->grid == VTGRID)) || (F->grid != ETAGRID) ) {
    field_error_message("Bad grid types to field2D_gradient_eta_to_uv()");
  }

  if (N != NFY) {
    field_error_message("Bad nx field2D sizes (F & FY) to field2D_gradient_eta_to_uv()");
  }

  if ((N-1) != NFX) {
    field_error_message("Bad nx field2D sizes (F & FX) to field2D_gradient_eta_to_uv()");
  }

  if ( (FY->M != M) && (FX->M != M)) {
    field_error_message("Bad M field2D sizes to field2D_gradient_eta_to_uv()");
  }


  for (i=0;i<N-1;i++) {
    fup = *(fd+M);
    fminus = *fd++;
    for (j=0;j<M-1;j++) {
      fplus = *fd++;
      *fyd++ = idy * (fplus - fminus);
      *fxd ++ = idx * (fup - fminus);
      fminus = fplus;
      fup = *(fd+MN);
    }
    *fxd++ = idx * (fup - fminus);
    *fyd++ = idy * ( *(fd-M) - fminus);
  }

  /* i = N-1 */

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


}

/* Calcualtes: FYY = \delta_y \delta_y (F) for either a VGRID OR ETAGRID */ 
/* Warning: Assumes that FYY is 2 X points smaller than F */

void field2D_diffy2(field2D *FYY, const field2D *F, double idy)
{
  int i,j;
  const int M = FYY->M;
  const int MN = M-1;
  const int NFYY = FYY->N;
  const int NF = F->N;

  const double *fd = &(F->data[M]);
  double *fyyd = &(FYY->data[0]);

  const double idy2 = idy*idy;

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

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

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

    /* j=0 */

    *fyyd++ = idy2 * ( *(fd+1) - 2.0 * (*fd) + *(fd+MN) );
    fd++;

    for (j=1;j<M-1;j++)  {
      *fyyd++ = idy2 * ( *(fd+1) - 2.0 * (*fd) + *(fd-1) );
      fd++;
    }

    /* j=M-1 */

    *fyyd++ = idy2 * ( *(fd-MN) - 2.0 * (*fd) + *(fd-1) );
    fd++;

  }
}


void field2D_laplacian(field2D *LAP_F, const field2D *F, const double idx, const double idy)
{

  int i,j;
  const int M = LAP_F->M;
  const int MN = M-1;
  const int NLAP_F = LAP_F->N;
  const int NF = F->N;

  const double *fd = &(F->data[M]);
  double *lapfd = &(LAP_F->data[0]);

  const double idy2 = idy*idy;

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

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

  field2D_diffx2(LAP_F, F, idx); /* LAP_F = \delta_x \delta_x (F) */

  /* Now add on the \delta_y \delta_y() part */

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

    /* j=0 */

    *lapfd++ += idy2 * ( *(fd+1) - 2.0 * (*fd) + *(fd+MN) );
    fd++;

    for (j=1;j<M-1;j++)  {
      *lapfd++ += idy2 * ( *(fd+1) - 2.0 * (*fd) + *(fd-1) );
      fd++;
    }

    /* j=M-1 */

    *lapfd++ += idy2 * ( *(fd-MN) - 2.0 * (*fd) + *(fd-1) );
    fd++;

  }
}






void field2D_apply_nogradient_bc(field2D *T)
{
  int N = T->N;
  int M = T->M;
  int j;

  double *td0 = &(T->data[0]);
  double *tdN = &(T->data[(N-1)*M]);

  for (j=0;j<M;j++) {
    td0[j] = td0[j+M];   // shore boundary condition
    tdN[j] = tdN[j-M];   // far slip BC
  }

}


void field2D_apply_zero_bc(field2D *T)
{
  int N = T->N;
  int M = T->M;
  int j;

  double *td0 = &(T->data[0]);
  double *tdN = &(T->data[(N-1)*M]);

  for (j=0;j<M;j++) {
    td0[j] = 0.0;     // shore boundary condition
    tdN[j] = 0.0;     // far slip BC
  }

}


/* ------------------------------------------ 
     Threads functions 
   ------------------------------------------ */

typedef struct {
  field2D *T;   /* the target */
  field2D *A;   /* the source  */
  field2D *B;   /* the 2nd source if needed */
  double idx;
  double idy;
  int istartT;
  int iendT;
  int istartA;
  int iendA;
} field2D_threads_t;


field2D_threads_t *D1;
field2D_threads_t *D2;  /* later can make this an pointer to array of D pointers for > 2 cpus */


pthread_t field2D_main_thread;
pthread_t work_thread1;
pthread_t work_thread2;  /* same here make it multiple work threads for > 2 cpus */

void *work_thread1_result;
void *work_thread2_result;
int num_cpu;



int field2D_copy_threads_work(field2D_threads_t *D)           /* T = A */
{

  field2D *T = D->T;
  field2D *A = D->A;

  int istartT = D->istartT;
  int istartA = D->istartA;
  int iendT = D->iendT;

  int M = T->M;
  int NUM = (iendT - istartT)*M;
  int i;

  double *td = &(T->data[istartT*M]);
  const double *ad = &(A->data[istartA*M]);

  
  /* Then go through and do the copy */

  for (i=0;i<NUM;i++)
    *td++ = *ad++;

  // pthread_exit(NULL);

}


int field2D_self_add_threads_work(field2D_threads_t *D)           /* T = A */
{

  field2D *T = D->T;
  field2D *A = D->A;

  int istartT = D->istartT;
  int istartA = D->istartA;
  int iendT = D->iendT;

  int M = T->M;
  int NUM = (iendT - istartT)*M;
  int i;

  double *td = &(T->data[istartT*M]);
  const double *ad = &(A->data[istartA*M]);

  //  fprintf(stderr,"** Starting field2D_self_add_threads_work() w/ istart=%d\n",istartT);

  
  /* Then go through and do the self add */

  for (i=0;i<NUM;i++)
    *td++ += *ad++;

  //  pthread_exit(NULL);

}



void field2D_copy_threads(field2D *T, const field2D *A)           /* T = A */
{

  int N = T->N;
  int M = T->M;

  gint status1, status2;  

  int ihalf;

  /* Check that the dimensions are ok */

  if ( A->M != M  ) {
    field_error_message("Bad field2D M sizes to void field2D_copy_threads()");
  }

  if ( !( ( A->N == N) || (A->N == N-2) ) ) {
    field_error_message("Bad field2D N sizes to void field2D_copy_threads()");
  }


  D1->T = T;
  D1->A = A;

  D2->T = T;
  D2->A = A;

  switch (A->N - N) {
  case 0:           /* if A & T have same sizes, OK, start at same location */
    ihalf = N/2;
    D1->istartA = 0;
    D1->istartT = 0;
    D1->iendA = ihalf;
    D1->iendT = ihalf;
    D2->istartA = ihalf;
    D2->istartT = ihalf;
    D2->iendA = T->N;
    D2->iendT = T->N;
    break;
  case -2:             /* if A(N,M) & T(N+2,M) then move td pointer up */
    ihalf = (A->N)/2;
    D1->istartA = 0;
    D1->istartT = 1;
    D1->iendA = ihalf;
    D1->iendT = ihalf+1;
    D2->istartA = ihalf;
    D2->istartT = ihalf+1;
    D2->iendA = A->N;
    D2->iendT = T->N-1;
    break;
  default:  field_error_message("Bad case to switch() {...}  in field2D_copy()");
  }



  status1 = pthread_create(&work_thread1, NULL, (void *) field2D_copy_threads_work, D1);
  if (status1 != 0)
    field_error_message("** Couldn't create thread!");

  status2 = pthread_create(&work_thread2, NULL, (void *) field2D_copy_threads_work, D2);
  if (status2 != 0)
    field_error_message("** Couldn't create thread!");

  pthread_join(work_thread1, &work_thread1_result);
  pthread_join(work_thread2, &work_thread2_result);


}


void field2D_self_add_threads(field2D *T, const field2D *A)           /* T += A */
{

  int N = T->N;
  int M = T->M;

  gint status1, status2;  

  int ihalf;

  /* Check that the dimensions are ok */

  if ( A->M != M  ) {
    field_error_message("Bad field2D M sizes to void field2D_self_add_threads()");
  }

  if ( !( ( A->N == N) || (A->N == N-2) ) ) {
    field_error_message("Bad field2D N sizes to void field2D_self_add_threads()");
  }


  D1->T = T;
  D1->A = A;

  D2->T = T;
  D2->A = A;

  switch (A->N - N) {
  case 0:           /* if A & T have same sizes, OK, start at same location */
    ihalf = N/2;
    D1->istartA = 0;
    D1->istartT = 0;
    D1->iendA = ihalf;
    D1->iendT = ihalf;
    D2->istartA = ihalf;
    D2->istartT = ihalf;
    D2->iendA = T->N;
    D2->iendT = T->N;
    break;
  case -2:             /* if A(N,M) & T(N+2,M) then move td pointer up */
    ihalf = (A->N)/2;
    D1->istartA = 0;
    D1->istartT = 1;
    D1->iendA = ihalf;
    D1->iendT = ihalf+1;
    D2->istartA = ihalf;
    D2->istartT = ihalf+1;
    D2->iendA = A->N;
    D2->iendT = T->N-1;
    break;
  default:  field_error_message("Bad case to switch() {...}  in field2D_self_add_threads()");
  }

  //  fprintf(stderr,"ihalf = %d\n",ihalf);


  status1 = pthread_create(&work_thread1, NULL, (void *) field2D_self_add_threads_work, D1);
  if (status1 != 0)
    field_error_message("** Couldn't create thread!");

  status2 = pthread_create(&work_thread2, NULL, (void *) field2D_self_add_threads_work, D2);
  if (status2 != 0)
    field_error_message("** Couldn't create thread!");

  pthread_join(work_thread1, NULL);
  pthread_join(work_thread2, NULL);


}

void field2D_init_threads(int numcpu)
{
  num_cpu = numcpu;
  D1 = (field2D_threads_t *) g_malloc(sizeof(field2D_threads_t));
  D2 = (field2D_threads_t *) g_malloc(sizeof(field2D_threads_t));


  field2D_main_thread = pthread_self();



}


void field2D_close_threads();