WIP: mpas particles

hguo · Oct 16, 2023 · bb7fb5c · bb7fb5c
1 parent ea38c76
commit bb7fb5c
Show file tree

Hide file tree

Showing 2 changed files with 53 additions and 12 deletions.
diff --git a/include/ftk/filters/mpas_ocean_particle_tracker.cuh b/include/ftk/filters/mpas_ocean_particle_tracker.cuh
@@ -27,12 +27,13 @@ typedef struct {
 
   // particle data
   int nparticles;
-  ftk::feature_point_lite_t *hcps = NULL, *dcps = NULL;
+  ftk::feature_point_lite_t *htrajs = NULL, *dtrajs = NULL;
 
 } mop_ctx_t;
 
 void mop_create_ctx(mop_ctx_t **c_, int device=0);
 void mop_destroy_ctx(mop_ctx_t **c_);
+
 void mop_load_mesh(mop_ctx_t *c,
     const int ncells, 
     const int nlayers, 
@@ -49,9 +50,11 @@ void mop_load_data(mop_ctx_t *c,
     const double *zTop, 
     const double *A);
 
-void mop_load_particles(mop_ctx_t *c,
+void mop_load_particles(mop_ctx_t *c, 
+    const int n,
+    ftk::feature_point_lite_t *);
 
-void mop_execute(mop_ctx_t *c, int scope, int current_timestep);
+void mop_execute(mop_ctx_t *c, int current_timestep);
 void mop_swap(mop_ctx_t *c);
 
 #endif
diff --git a/src/filters/particle_tracer_mpas.cu b/src/filters/particle_tracer_mpas.cu
@@ -227,12 +227,50 @@ void mop_destroy_ctx(mop_ctx_t **c_)
   ctx_t *c = *c_;
 
   if (c->d_Xc != NULL) cudaFree(c->d_Xc);
-  // TODO
+  if (c->d_Xv != NULL) cudaFree(c->d_Xv);
+  if (c->d_nedges_on_cell != NULL) cudaFree(c->d_nedges_on_cell);
+  if (c->d_cells_on_cell != NULL) cudaFree(c->d_cells_on_cell);
+  if (c->d_verts_on_cell != NULL) cudaFree(c->d_verts_on_cell);
+
+  if (c->d_V[0] != NULL) cudaFree(c->d_V[0]);
+  if (c->d_V[1] != NULL) cudaFree(c->d_V[1]);
+  if (c->d_Vv[0] != NULL) cudaFree(c->d_Vv[0]);
+  if (c->d_Vv[1] != NULL) cudaFree(c->d_Vv[1]);
+  if (c->d_zTop[0] != NULL) cudaFree(c->d_zTop[0]);
+  if (c->d_zTop[1] != NULL) cudaFree(c->d_zTop[1]);
+  if (c->d_A[0] != NULL) cudaFree(c->d_A[0]);
+  if (c->d_A[1] != NULL) cudaFree(c->d_A[1]);
 
   free(*c_);
   *c_ = NULL;
 }
 
+void mop_load_particles(mop_ctx_t *c, 
+    const int n, 
+    ftk::feature_point_lite_t *data)
+{
+  if (c->htrajs == NULL) {
+    c->htrajs = (ftk::feature_point_lite_t*)
+      malloc(n * sizeof(ftk::feature_point_lite_t));
+  } else {
+    c->htrajs = (ftk::feature_point_lite_t*)
+      realloc(c->htrajs, n * sizeof(ftk::feature_point_lite_t));
+  }
+
+  if (c->dtrajs == NULL) {
+    cudaMalloc((void**)&c->dtrajs, 
+        n * sizeof(ftk::feature_point_lite_t));
+  } else {
+    if (n != c->nparticles) {
+      cudaFree(c->dtrajs);
+      cudaMalloc((void**)&c->dtrajs, 
+          n * sizeof(ftk::feature_point_lite_t));
+    }
+  }
+
+  c->nparticles = n;
+}
+
 void mop_load_mesh(mop_ctx_t *c,
     const int ncells, 
     const int nlayers, 
@@ -306,12 +344,9 @@ void mop_load_data(mop_ctx_t *c,
 }
 
 
-void mop_execute(mop_ctx_t *c, int scope, int current_timestep)
+void mop_execute(mop_ctx_t *c, int current_timestep)
 {
-  size_t ntasks;
-  if (scope == scope_ordinal) ntasks = mx4n2_ordinal;
-  else ntasks = mx4n2_interval;
-
+  size_t ntasks = c->nparticles;
   fprintf(stderr, "ntasks=%zu\n", ntasks);
 
   const int maxGridDim = 1024;
@@ -322,6 +357,8 @@ void mop_execute(mop_ctx_t *c, int scope, int current_timestep)
   if (nBlocks >= maxGridDim) gridSize = dim3(idivup(nBlocks, maxGridDim), maxGridDim);
   else gridSize = dim3(nBlocks);
 
+  // TODO
+#if 0
   sweep_simplices<int, double><<<gridSize, blockSize>>>(
       scope, current_timestep, 
       c->factor,
@@ -344,12 +381,13 @@ void mop_execute(mop_ctx_t *c, int scope, int current_timestep)
   cudaMemcpy(c->hcps, c->dcps, sizeof(cp_t) * c->hncps, cudaMemcpyDeviceToHost);
 
   checkLastCudaError("[FTK-CUDA] cuda memcpy device to host, 2");
+#endif
 }
 
 void mop_swap(mop_ctx_t *c)
 {
   std::swap(c->d_V[0], c->d_V[1]);
-  std::swap(c->d_vector[0], c->d_vector[1]);
-  std::swap(c->d_jacobian[0], c->d_jacobian[1]);
+  std::swap(c->d_Vv[0], c->d_Vv[1]);
+  std::swap(c->d_zTop[0], c->d_zTop[1]);
+  std::swap(c->d_A[0], c->d_A[1]);
 }
-#endif