cudaThreadExit() was not being called, which caused failures on some platforms

platforms/cuda/src/kernels/gpu.cpp

@@ -1692,11 +1692,8 @@ void* gpuInit(int numAtoms, unsigned int device, bool useBlockingSync)
     int SMMinor = 0;
 
     // Select which device to use
-    int currentDevice;
-    cudaError_t status = cudaGetDevice(&currentDevice);
-    RTERROR(status, "Error getting CUDA device")
-    if (device != currentDevice)
-        cudaSetDevice(device); // Ignore errors
+    cudaError_t status = cudaSetDevice(device);
+    RTERROR(status, "Error setting CUDA device")
     status = cudaGetDevice(&gpu->device);
     RTERROR(status, "Error getting CUDA device")
     status = cudaSetDeviceFlags(useBlockingSync ? cudaDeviceBlockingSync : cudaDeviceScheduleAuto);
@@ -1705,7 +1702,7 @@ void* gpuInit(int numAtoms, unsigned int device, bool useBlockingSync)
 
     // Determine kernel call configuration
     cudaDeviceProp deviceProp;
-    cudaGetDeviceProperties(&deviceProp, currentDevice);
+    cudaGetDeviceProperties(&deviceProp, gpu->device);
 
     // Determine SM version
     if (deviceProp.major == 1)
@@ -2089,6 +2086,7 @@ void gpuShutDown(gpuContext gpu)
 
     // Wrap up
     delete gpu;
+    cudaThreadExit();
     return;
 }
 

platforms/cuda/tests/TestCudaCustomNonbondedForce.cpp

@@ -294,17 +294,27 @@ void testCoulombLennardJones() {
     customSystem.addForce(customNonbonded);
     VerletIntegrator integrator1(0.01);
     VerletIntegrator integrator2(0.01);
-    Context context1(standardSystem, integrator1, platform);
-    context1.setPositions(positions);
-    context1.setVelocities(velocities);
-    State state1 = context1.getState(State::Forces | State::Energy);
-    Context context2(customSystem, integrator2, platform);
-    context2.setPositions(positions);
-    context2.setVelocities(velocities);
-    State state2 = context2.getState(State::Forces | State::Energy);
-    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-4);
+    double energy1, energy2;
+    vector<Vec3> forces1, forces2;
+    {
+        Context context(standardSystem, integrator1, platform);
+        context.setPositions(positions);
+        context.setVelocities(velocities);
+        State state = context.getState(State::Forces | State::Energy);
+        energy1 = state.getPotentialEnergy();
+        forces1 = state.getForces();
+    }
+    {
+        Context context(customSystem, integrator2, platform);
+        context.setPositions(positions);
+        context.setVelocities(velocities);
+        State state = context.getState(State::Forces | State::Energy);
+        energy2 = state.getPotentialEnergy();
+        forces2 = state.getForces();
+    }
+    ASSERT_EQUAL_TOL(energy1, energy2, 1e-4);
     for (int i = 0; i < numParticles; i++) {
-        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
+        ASSERT_EQUAL_VEC(forces1[i], forces2[i], 1e-4);
     }
 }
 

platforms/cuda/tests/TestCudaEwald.cpp

@@ -116,11 +116,11 @@ void testEwaldPME() {
         Vec3 f = cudaState.getForces()[i];
         positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
     }
-    Context cudaContext2(system, integrator, cuda);
-    cudaContext2.setPositions(positions);
+    cudaContext.reinitialize();
+    cudaContext.setPositions(positions);
     
     tol = 1e-3;
-    State cudaState2 = cudaContext2.getState(State::Energy);
+    State cudaState2 = cudaContext.getState(State::Energy);
     ASSERT_EQUAL_TOL(norm, (cudaState2.getPotentialEnergy()-cudaState.getPotentialEnergy())/delta, tol)
 
 //    (3)  Check whether the Reference and Cuda platforms agree when using PME
@@ -154,11 +154,11 @@ void testEwaldPME() {
         Vec3 f = cudaState.getForces()[i];
         positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
     }
-    Context cudaContext3(system, integrator, cuda);
-    cudaContext3.setPositions(positions);
+    cudaContext.reinitialize();
+    cudaContext.setPositions(positions);
      
     tol = 1e-3;
-    State cudaState3 = cudaContext3.getState(State::Energy);
+    State cudaState3 = cudaContext.getState(State::Energy);
     ASSERT_EQUAL_TOL(norm, (cudaState3.getPotentialEnergy()-cudaState.getPotentialEnergy())/delta, tol)
 }
 

platforms/cuda/tests/TestCudaNonbondedForce.cpp

@@ -131,7 +131,7 @@ void testExclusionsAnd14() {
     for (int i = 1; i < 5; ++i) {
  
         // Test LJ forces
-        
+
         vector<Vec3> positions(5);
         const double r = 1.0;
         for (int j = 0; j < 5; ++j) {
@@ -143,49 +143,55 @@ void testExclusionsAnd14() {
         nonbonded->setExceptionParameters(first14, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
         nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0.0);
         positions[i] = Vec3(r, 0, 0);
-        Context context(system, integrator, platform);
-        context.setPositions(positions);
-        State state = context.getState(State::Forces | State::Energy);
-        const vector<Vec3>& forces = state.getForces();
-        double x = 1.5/r;
-        double eps = 1.0;
-        double force = 4.0*eps*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/r;
-        double energy = 4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0));
-        if (i == 3) {
-            force *= 0.5;
-            energy *= 0.5;
-        }
-        if (i < 3) {
-            force = 0;
-            energy = 0;
+        // The following is in its own block, because CUDA can't deal with multiple Contexts
+        // existing on the same thread at the same time.
+        {
+            Context context(system, integrator, platform);
+            context.setPositions(positions);
+            State state = context.getState(State::Forces | State::Energy);
+            const vector<Vec3>& forces = state.getForces();
+            double x = 1.5/r;
+            double eps = 1.0;
+            double force = 4.0*eps*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/r;
+            double energy = 4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0));
+            if (i == 3) {
+                force *= 0.5;
+                energy *= 0.5;
+            }
+            if (i < 3) {
+                force = 0;
+                energy = 0;
+            }
+            ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
+            ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[i], TOL);
+            ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
         }
-        ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
-        ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[i], TOL);
-        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
 
         // Test Coulomb forces
-        
-        nonbonded->setParticleParameters(0, 2, 1.5, 0);
-        nonbonded->setParticleParameters(i, 2, 1.5, 0);
-        nonbonded->setExceptionParameters(first14, 0, 3, i == 3 ? 4/1.2 : 0, 1.5, 0);
-        nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0);
-        Context context2(system, integrator, platform);
-        context2.setPositions(positions);
-        state = context2.getState(State::Forces | State::Energy);
-        const vector<Vec3>& forces2 = state.getForces();
-        force = ONE_4PI_EPS0*4/(r*r);
-        energy = ONE_4PI_EPS0*4/r;
-        if (i == 3) {
-            force /= 1.2;
-            energy /= 1.2;
-        }
-        if (i < 3) {
-            force = 0;
-            energy = 0;
+
+        {
+            nonbonded->setParticleParameters(0, 2, 1.5, 0);
+            nonbonded->setParticleParameters(i, 2, 1.5, 0);
+            nonbonded->setExceptionParameters(first14, 0, 3, i == 3 ? 4/1.2 : 0, 1.5, 0);
+            nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0);
+            Context context(system, integrator, platform);
+            context.setPositions(positions);
+            State state = context.getState(State::Forces | State::Energy);
+            const vector<Vec3>& forces2 = state.getForces();
+            double force = ONE_4PI_EPS0*4/(r*r);
+            double energy = ONE_4PI_EPS0*4/r;
+            if (i == 3) {
+                force /= 1.2;
+                energy /= 1.2;
+            }
+            if (i < 3) {
+                force = 0;
+                energy = 0;
+            }
+            ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces2[0], TOL);
+            ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces2[i], TOL);
+            ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
         }
-        ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces2[0], TOL);
-        ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces2[i], TOL);
-        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
     }
 }