diff --git a/src/CBDisk.cc b/src/CBDisk.cc index ee8d9836..67137cab 100644 --- 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 index 9d512b32..520611a7 100644 --- 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::max(), bool last=false); + void Run(int nstep, int mstep=std::numeric_limits::max(), bool final=false); }; #endif diff --git a/src/PolarBasis.cc b/src/PolarBasis.cc index 89beb1fe..aae196ae 100644 --- 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 index 468a7888..5a579020 100644 --- a/src/cudaBiorthCyl.cu +++ b/src/cudaBiorthCyl.cu @@ -187,14 +187,19 @@ void BiorthCyl::initialize_cuda // Add background arrays // std::vector> 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; ipot(r); - tt[1][i] = disk->dpot(r); + for (int i=0; ipot(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 index 2f7c0afb..008a03e7 100644 --- 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);