<commit_before>#include <stdio.h>
#include <mpi.h>
#include "fem.h"

#include <occa.hpp>

occa::device device;

occa::memory c_LIFT;
occa::memory c_DrDsDt;
occa::memory c_surfinfo;
occa::memory c_mapinfo;
occa::memory c_vgeo;
occa::memory c_Q; 
occa::memory c_partQ; 
occa::memory c_rhsQ; 
occa::memory c_resQ; 
occa::memory c_tmp;
occa::memory c_parmapOUT;

occa::kernel rkKernel;
occa::kernel volumeKernel;
occa::kernel surfaceKernel;
occa::kernel partialGetKernel;

template <typename T>
void diagnose_array(const char *message, occa::memory &c_a, int N){

  device.finish();
  
  T *h_a = (T*) calloc(N, sizeof(T));
  
  c_a.copyTo(h_a, N*sizeof(T));

  double suma = 0;
  for(int n=0;n<N;++n){
    suma += h_a[n];
  }

  printf("%s: sum = %17.15g\n", message, suma);
  
  free(h_a);
}




double InitOCCA3d(Mesh *mesh, int Nfields){
  
  int plat = 1, dev = 0;
  device.setup("OpenCL", plat, dev);

  /* Q  */
  int sz = mesh->K*(BSIZE)*p_Nfields*sizeof(float); 

  float *f_Q = (float*) calloc(mesh->K*BSIZE*p_Nfields, sizeof(float));

  c_Q    = device.malloc(sz, f_Q);
  c_rhsQ = device.malloc(sz, f_Q);
  c_resQ = device.malloc(sz, f_Q);

  printf("sz1= %d\n", sz);
  
  sz = mesh->parNtotalout*sizeof(float);
  c_tmp  = device.malloc(sz, f_Q); // should not use f_Q
  c_partQ = device.malloc(sz, f_Q);

  printf("sz2= %d\n", sz);

  /*  LIFT  */
  sz = p_Np*(p_Nfp)*p_Nfaces*sizeof(float);

  float *f_LIFT = (float*) malloc(sz);
  int skL = 0;
  for(int m=0;m<p_Nfp;++m){
    for(int n=0;n<p_Np;++n){
      for(int f=0;f<p_Nfaces;++f){
	f_LIFT[skL++] = mesh->LIFT[0][p_Nfp*p_Nfaces*n+(f+p_Nfaces*m)];
      }
    }
  }

  c_LIFT = device.malloc(sz, f_LIFT);
   
  /* DrDsDt */
  sz = BSIZE*BSIZE*4*sizeof(float);

  float* h_DrDsDt = (float*) calloc(BSIZE*BSIZE*4, sizeof(float));
  int sk = 0;
  /* note transposed arrays to avoid "bank conflicts" */
  for(int n=0;n<p_Np;++n){
    for(int m=0;m<p_Np;++m){
      h_DrDsDt[4*(m+n*BSIZE)+0] = mesh->Dr[0][n+m*p_Np];
      h_DrDsDt[4*(m+n*BSIZE)+1] = mesh->Ds[0][n+m*p_Np];
      h_DrDsDt[4*(m+n*BSIZE)+2] = mesh->Dt[0][n+m*p_Np];
    }
  }
   
  c_DrDsDt = device.malloc(sz, h_DrDsDt);
   
  free(h_DrDsDt);

  /* vgeo */
  double drdx, dsdx, dtdx;
  double drdy, dsdy, dtdy;
  double drdz, dsdz, dtdz, J;
  float *vgeo = (float*) calloc(12*mesh->K, sizeof(float));

  for(int k=0;k<mesh->K;++k){
    GeometricFactors3d(mesh, k, 
		       &drdx, &dsdx, &dtdx,
		       &drdy, &dsdy, &dtdy,
		       &drdz, &dsdz, &dtdz, &J);

    vgeo[k*12+0] = drdx; vgeo[k*12+1] = drdy; vgeo[k*12+2] = drdz;
    vgeo[k*12+4] = dsdx; vgeo[k*12+5] = dsdy; vgeo[k*12+6] = dsdz;
    vgeo[k*12+8] = dtdx; vgeo[k*12+9] = dtdy; vgeo[k*12+10] = dtdz;

  }

  sz = mesh->K*12*sizeof(float);
  c_vgeo = device.malloc(sz, vgeo);
   
  /* surfinfo (vmapM, vmapP, Fscale, Bscale, nx, ny, nz, 0) */
  int sz5 = mesh->K*p_Nfp*p_Nfaces*5*sizeof(float); 
  float* h_surfinfo = (float*) malloc(sz5); 

  int sz2 = mesh->K*p_Nfp*p_Nfaces*2*sizeof(int); 
  int* h_mapinfo = (int*) malloc(sz2); 
   
  /* local-local info */
  sk = 0;
  int skP = -1;
  double *nxk = BuildVector(mesh->Nfaces);
  double *nyk = BuildVector(mesh->Nfaces);
  double *nzk = BuildVector(mesh->Nfaces);
  double *sJk = BuildVector(mesh->Nfaces);

  double dt = 1e6;

  for(int k=0;k<mesh->K;++k){

    GeometricFactors3d(mesh, k, 
		       &drdx, &dsdx, &dtdx,
		       &drdy, &dsdy, &dtdy,
		       &drdz, &dsdz, &dtdz, &J);

    Normals3d(mesh, k, nxk, nyk, nzk, sJk);
     
    for(int f=0;f<mesh->Nfaces;++f){

      dt = min(dt, J/sJk[f]);
       
      for(int m=0;m<p_Nfp;++m){
	int n = m + f*p_Nfp + p_Nfp*p_Nfaces*k;
	int idM = mesh->vmapM[n];
	int idP = mesh->vmapP[n];
	int  nM = idM%p_Np; 
	int  nP = idP%p_Np; 
	int  kM = (idM-nM)/p_Np;
	int  kP = (idP-nP)/p_Np;
	idM = nM + Nfields*BSIZE*kM;
	idP = nP + Nfields*BSIZE*kP;
	 
	/* stub resolve some other way */
	if(mesh->vmapP[n]<0){
	  idP = mesh->vmapP[n]; /* -ve numbers */
	}
 
	sk = 2*p_Nfp*p_Nfaces*k+m+f*p_Nfp;
	h_mapinfo[sk + 0*p_Nfp*p_Nfaces] = idM;
	h_mapinfo[sk + 1*p_Nfp*p_Nfaces] = idP;

	sk = 5*p_Nfp*p_Nfaces*k+m+f*p_Nfp;
	h_surfinfo[sk + 0*p_Nfp*p_Nfaces] = sJk[f]/(2.*J);
	h_surfinfo[sk + 1*p_Nfp*p_Nfaces] = (idM==idP)?-1.:1.;
	h_surfinfo[sk + 2*p_Nfp*p_Nfaces] = nxk[f];
	h_surfinfo[sk + 3*p_Nfp*p_Nfaces] = nyk[f];
	h_surfinfo[sk + 4*p_Nfp*p_Nfaces] = nzk[f];
      }
    }
  }
   
  c_mapinfo = device.malloc(sz2, h_mapinfo);
  c_surfinfo = device.malloc(sz5, h_surfinfo);

  free(h_mapinfo);
  free(h_surfinfo);

  printf("mesh->parNtotalout=%d\n", mesh->parNtotalout);
  sz = mesh->parNtotalout*sizeof(int);
  c_parmapOUT = device.malloc(sz, mesh->parmapOUT);

  /* now build kernels */
  occa::kernelInfo dgInfo;
   
  dgInfo.addDefine("p_Np",      p_Np);
  dgInfo.addDefine("p_Nfp",     p_Nfp);
  dgInfo.addDefine("p_Nfaces",  p_Nfaces);
  dgInfo.addDefine("p_Nfields", p_Nfields);
  dgInfo.addDefine("BSIZE",     BSIZE);

  int Nt, Nb;
  {
    volumeKernel = device.buildKernelFromSource("src/MaxwellsVolumeKernel3D.occa", 
						"MaxwellsVolumeKernel3D",
						dgInfo);
    Nt = p_Np;
    Nb = mesh->K;
    volumeKernel.setWorkingDims(1, Nt, Nb);
  }
  
  {<commit_msg> padded parallel buffers by 1<commit_after>#include <stdio.h>
#include <mpi.h>
#include "fem.h"

#include <occa.hpp>

occa::device device;

occa::memory c_LIFT;
occa::memory c_DrDsDt;
occa::memory c_surfinfo;
occa::memory c_mapinfo;
occa::memory c_vgeo;
occa::memory c_Q; 
occa::memory c_partQ; 
occa::memory c_rhsQ; 
occa::memory c_resQ; 
occa::memory c_tmp;
occa::memory c_parmapOUT;

occa::kernel rkKernel;
occa::kernel volumeKernel;
occa::kernel surfaceKernel;
occa::kernel partialGetKernel;

template <typename T>
void diagnose_array(const char *message, occa::memory &c_a, int N){

  device.finish();
  
  T *h_a = (T*) calloc(N, sizeof(T));
  
  c_a.copyTo(h_a, N*sizeof(T));

  double suma = 0;
  for(int n=0;n<N;++n){
    suma += h_a[n];
  }

  printf("%s: sum = %17.15g\n", message, suma);
  
  free(h_a);
}




double InitOCCA3d(Mesh *mesh, int Nfields){
  
  int plat = 0, dev = 2;
  device.setup("OpenCL", plat, dev);

  /* Q  */
  int sz = mesh->K*(BSIZE)*p_Nfields*sizeof(float); 

  float *f_Q = (float*) calloc(mesh->K*BSIZE*p_Nfields, sizeof(float));

  c_Q    = device.malloc(sz, f_Q);
  c_rhsQ = device.malloc(sz, f_Q);
  c_resQ = device.malloc(sz, f_Q);

  printf("sz1= %d\n", sz);
  
  sz = mesh->parNtotalout*sizeof(float);
  c_tmp  = device.malloc(sz+1, f_Q); // should not use f_Q
  c_partQ = device.malloc(sz+1, f_Q);

  printf("sz2= %d\n", sz);

  /*  LIFT  */
  sz = p_Np*(p_Nfp)*p_Nfaces*sizeof(float);

  float *f_LIFT = (float*) malloc(sz);
  int skL = 0;
  for(int m=0;m<p_Nfp;++m){
    for(int n=0;n<p_Np;++n){
      for(int f=0;f<p_Nfaces;++f){
	f_LIFT[skL++] = mesh->LIFT[0][p_Nfp*p_Nfaces*n+(f+p_Nfaces*m)];
      }
    }
  }

  c_LIFT = device.malloc(sz, f_LIFT);
   
  /* DrDsDt */
  sz = BSIZE*BSIZE*4*sizeof(float);

  float* h_DrDsDt = (float*) calloc(BSIZE*BSIZE*4, sizeof(float));
  int sk = 0;
  /* note transposed arrays to avoid "bank conflicts" */
  for(int n=0;n<p_Np;++n){
    for(int m=0;m<p_Np;++m){
      h_DrDsDt[4*(m+n*BSIZE)+0] = mesh->Dr[0][n+m*p_Np];
      h_DrDsDt[4*(m+n*BSIZE)+1] = mesh->Ds[0][n+m*p_Np];
      h_DrDsDt[4*(m+n*BSIZE)+2] = mesh->Dt[0][n+m*p_Np];
    }
  }
   
  c_DrDsDt = device.malloc(sz, h_DrDsDt);
   
  free(h_DrDsDt);

  /* vgeo */
  double drdx, dsdx, dtdx;
  double drdy, dsdy, dtdy;
  double drdz, dsdz, dtdz, J;
  float *vgeo = (float*) calloc(12*mesh->K, sizeof(float));

  for(int k=0;k<mesh->K;++k){
    GeometricFactors3d(mesh, k, 
		       &drdx, &dsdx, &dtdx,
		       &drdy, &dsdy, &dtdy,
		       &drdz, &dsdz, &dtdz, &J);

    vgeo[k*12+0] = drdx; vgeo[k*12+1] = drdy; vgeo[k*12+2] = drdz;
    vgeo[k*12+4] = dsdx; vgeo[k*12+5] = dsdy; vgeo[k*12+6] = dsdz;
    vgeo[k*12+8] = dtdx; vgeo[k*12+9] = dtdy; vgeo[k*12+10] = dtdz;

  }

  sz = mesh->K*12*sizeof(float);
  c_vgeo = device.malloc(sz, vgeo);
   
  /* surfinfo (vmapM, vmapP, Fscale, Bscale, nx, ny, nz, 0) */
  int sz5 = mesh->K*p_Nfp*p_Nfaces*5*sizeof(float); 
  float* h_surfinfo = (float*) malloc(sz5); 

  int sz2 = mesh->K*p_Nfp*p_Nfaces*2*sizeof(int); 
  int* h_mapinfo = (int*) malloc(sz2); 
   
  /* local-local info */
  sk = 0;
  int skP = -1;
  double *nxk = BuildVector(mesh->Nfaces);
  double *nyk = BuildVector(mesh->Nfaces);
  double *nzk = BuildVector(mesh->Nfaces);
  double *sJk = BuildVector(mesh->Nfaces);

  double dt = 1e6;

  for(int k=0;k<mesh->K;++k){

    GeometricFactors3d(mesh, k, 
		       &drdx, &dsdx, &dtdx,
		       &drdy, &dsdy, &dtdy,
		       &drdz, &dsdz, &dtdz, &J);

    Normals3d(mesh, k, nxk, nyk, nzk, sJk);
     
    for(int f=0;f<mesh->Nfaces;++f){

      dt = min(dt, J/sJk[f]);
       
      for(int m=0;m<p_Nfp;++m){
	int n = m + f*p_Nfp + p_Nfp*p_Nfaces*k;
	int idM = mesh->vmapM[n];
	int idP = mesh->vmapP[n];
	int  nM = idM%p_Np; 
	int  nP = idP%p_Np; 
	int  kM = (idM-nM)/p_Np;
	int  kP = (idP-nP)/p_Np;
	idM = nM + Nfields*BSIZE*kM;
	idP = nP + Nfields*BSIZE*kP;
	 
	/* stub resolve some other way */
	if(mesh->vmapP[n]<0){
	  idP = mesh->vmapP[n]; /* -ve numbers */
	}
 
	sk = 2*p_Nfp*p_Nfaces*k+m+f*p_Nfp;
	h_mapinfo[sk + 0*p_Nfp*p_Nfaces] = idM;
	h_mapinfo[sk + 1*p_Nfp*p_Nfaces] = idP;

	sk = 5*p_Nfp*p_Nfaces*k+m+f*p_Nfp;
	h_surfinfo[sk + 0*p_Nfp*p_Nfaces] = sJk[f]/(2.*J);
	h_surfinfo[sk + 1*p_Nfp*p_Nfaces] = (idM==idP)?-1.:1.;
	h_surfinfo[sk + 2*p_Nfp*p_Nfaces] = nxk[f];
	h_surfinfo[sk + 3*p_Nfp*p_Nfaces] = nyk[f];
	h_surfinfo[sk + 4*p_Nfp*p_Nfaces] = nzk[f];
      }
    }
  }
   
  c_mapinfo = device.malloc(sz2, h_mapinfo);
  c_surfinfo = device.malloc(sz5, h_surfinfo);

  free(h_mapinfo);
  free(h_surfinfo);

  printf("mesh->parNtotalout=%d\n", mesh->parNtotalout);
  sz = mesh->parNtotalout*sizeof(int);
  c_parmapOUT = device.malloc(sz+1, mesh->parmapOUT);

  /* now build kernels */
  occa::kernelInfo dgInfo;
   
  dgInfo.addDefine("p_Np",      p_Np);
  dgInfo.addDefine("p_Nfp",     p_Nfp);
  dgInfo.addDefine("p_Nfaces",  p_Nfaces);
  dgInfo.addDefine("p_Nfields", p_Nfields);
  dgInfo.addDefine("BSIZE",     BSIZE);

  int Nt, Nb;
  {
    volumeKernel = device.buildKernelFromSource("src/MaxwellsVolumeKernel3D.occa", 
						"MaxwellsVolumeKernel3D",
						dgInfo);
    Nt = p_Np;
    Nb = mesh->K;
    volumeKernel.setWorkingDims(1, Nt, Nb);
  }
  
  {
{'fraction_non_alphanumeric': 0.1942564300908915, 'fraction_numerical': 0.01933861922258751, 'mean_word_length': 3.556387665198238, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 65, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7112908', 'n_tokens_mistral': 5305, 'n_tokens_neox': 5035, 'n_words': 825}