Merge branch 'VTKfiles' of github.com:EXP-code/EXP into VTKfiles

EXP-code · Nov 21, 2024 · 34a8dcf · 34a8dcf
2 parents fac5609 + 6e5dbe4
commit 34a8dcf
Show file tree

Hide file tree

Showing 5 changed files with 57 additions and 15 deletions.
diff --git a/src/CBDisk.cc b/src/CBDisk.cc
@@ -20,6 +20,7 @@ CBDisk::CBDisk(Component* c0, const YAML::Node& conf, MixtureBasis* m) :
   PolarBasis(c0, conf, m)
 {
   id = "CBDisk";
+  is_flat = true;
 
   // Radial scale factor/softening radius
   //

diff --git a/src/OutputContainer.H b/src/OutputContainer.H
@@ -32,7 +32,7 @@ public:
   void initialize();
 
   //! Execute the all methods in the container
-  void Run(int nstep, int mstep=std::numeric_limits<int>::max(), bool last=false);
+  void Run(int nstep, int mstep=std::numeric_limits<int>::max(), bool final=false);
 };
 
 #endif

diff --git a/src/PolarBasis.cc b/src/PolarBasis.cc
@@ -401,8 +401,20 @@ void * PolarBasis::determine_coefficients_thread(void * arg)
 {
   // For biorthogonal density component and normalization
   // 
-  constexpr double norm0 = 1.0;
-  constexpr double norm1 = M_SQRT2;
+  constexpr double norm0_3d = 2.0*M_PI * 0.5*M_2_SQRTPI/M_SQRT2;
+  constexpr double norm1_3d = 2.0*M_PI * 0.5*M_2_SQRTPI;
+
+  constexpr double norm0_2d = 1.0;
+  constexpr double norm1_2d = M_SQRT2;
+
+  double norm0, norm1;
+  if (is_flat) {
+    norm0 = norm0_2d;
+    norm1 = norm1_2d;
+  } else {
+    norm0 = norm0_3d;
+    norm1 = norm1_3d;
+  }
 
   double r, r2, facL=1.0, fac1, fac2, phi, mass;
   double xx, yy, zz;
@@ -1044,8 +1056,20 @@ void PolarBasis::multistep_update(int from, int to, Component *c, int i, int id)
 
   // For biorthogonal density component and normalization
   // 
-  constexpr double norm0 = 1.0;
-  constexpr double norm1 = M_SQRT2;
+  constexpr double norm0_3d = 2.0*M_PI * 0.5*M_2_SQRTPI/M_SQRT2;
+  constexpr double norm1_3d = 2.0*M_PI * 0.5*M_2_SQRTPI;
+
+  constexpr double norm0_2d = 1.0;
+  constexpr double norm1_2d = M_SQRT2;
+
+  double norm0, norm1;
+  if (is_flat) {
+    norm0 = norm0_2d;
+    norm1 = norm1_2d;
+  } else {
+    norm0 = norm0_3d;
+    norm1 = norm1_3d;
+  }
 
   double mass = c->Mass(i) * component->Adiabatic();
 
@@ -1323,8 +1347,20 @@ void PolarBasis::compute_multistep_coefficients()
 
 void * PolarBasis::determine_acceleration_and_potential_thread(void * arg)
 {
-  constexpr double norm0 = 1.0;
-  constexpr double norm1 = M_SQRT2;
+  constexpr double norm0_3d = 0.5*M_2_SQRTPI/M_SQRT2;
+  constexpr double norm1_3d = 0.5*M_2_SQRTPI;
+
+  constexpr double norm0_2d = 1.0;
+  constexpr double norm1_2d = M_SQRT2;
+
+  double norm0, norm1;
+  if (is_flat) {
+    norm0 = norm0_2d;
+    norm1 = norm1_2d;
+  } else {
+    norm0 = norm0_3d;
+    norm1 = norm1_3d;
+  }
 
   double r, r0=0.0, phi;
   double potr, potz, potl, potp, p, pc, drc, drs, dzc, dzs, ps, dfacp, facdp;

diff --git a/src/cudaBiorthCyl.cu b/src/cudaBiorthCyl.cu
@@ -187,14 +187,19 @@ void BiorthCyl::initialize_cuda
   // Add background arrays
   //
   std::vector<thrust::host_vector<cuFP_t>> tt(ndim2);
-  for (auto & v : tt) v.resize(numr);
+  for (auto & v : tt) {
+    v.resize(numr);
+    thrust::fill(v.begin(), v.end(), 0.0);
+  }
 
-  double dx0  = (xmax - xmin)/(numr - 1);
+  if (disk) {
+    double dx0  = (xmax - xmin)/(numr - 1);
 
-  for (int i=0; i<numr; i++) {
-    double r = xi_to_r(xmin + dx0*i);
-    tt[0][i] = disk->pot(r);
-    tt[1][i] = disk->dpot(r);
+    for (int i=0; i<numr; i++) {
+      double r = xi_to_r(xmin + dx0*i);
+      tt[0][i] = disk->pot(r);
+      tt[1][i] = disk->dpot(r);
+    }
   }
 
   // Allocate CUDA array in device memory (a one-dimension 'channel')

diff --git a/src/end.cc b/src/end.cc
@@ -15,9 +15,9 @@
 void clean_up(void)
 {
 				// Call for final output to files
-  output->Run(this_step, 0, true);
+  if (output) output->Run(this_step, 0, true);
 				// Cache for restart
-  external->finish();
+  if (external) external->finish();
 
   MPI_Barrier(MPI_COMM_WORLD);