diff --git a/examples/thermal_deformation.cpp b/examples/thermal_deformation.cpp
index b008ad48..91b57163 100644
--- a/examples/thermal_deformation.cpp
+++ b/examples/thermal_deformation.cpp
@@ -98,7 +98,7 @@ int main( int argc, char* argv[] )
         particles->updateParticles( exec_space{}, init_functor );
 
         auto cabana_pd = CabanaPD::createSolverElastic<memory_space>(
-            inputs, particles, force_model );
+            inputs, particles, force_model, bc );
         cabana_pd->init_force();
         cabana_pd->run();
 
diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp
index c9973d8b..3136ca98 100644
--- a/src/CabanaPD_Boundary.hpp
+++ b/src/CabanaPD_Boundary.hpp
@@ -152,7 +152,7 @@ struct BoundaryCondition<BCIndexSpace, ZeroBCTag>
     }
 
     template <class ExecSpace, class ParticleType>
-    void apply( ExecSpace, ParticleType )
+    void apply( ExecSpace, ParticleType, double )
     {
     }
 };
@@ -188,8 +188,6 @@ struct BoundaryCondition<BCIndexSpace, ForceValueBCTag>
                 auto pid = index_space( b );
                 // This is specifically for the thermal deformation problem
                 temp( pid ) = 5000 * x( pid, 1 ) * t;
-                std::cout << "1: temp(" << pid << ")=" << temp( pid )
-                          << std::endl;
             } );
     }
 };
@@ -225,8 +223,6 @@ struct BoundaryCondition<BCIndexSpace, ForceUpdateBCTag>
                 auto pid = index_space( b );
                 // This is specifically for the thermal deformation problem
                 temp( pid ) += 5000 * x( pid, 1 ) * t;
-                std::cout << "2: temp(" << pid << ")=" << temp( pid )
-                          << std::endl;
             } );
     }
 };
@@ -263,8 +259,6 @@ struct BoundaryCondition<BCIndexSpace, ForceCrackBranchBCTag>
                 auto pid = index_space( b );
                 // This is specifically for the thermal deformation problem
                 temp( pid ) += 5000 * x( pid, 1 ) * t;
-                std::cout << "3: temp(" << pid << ")=" << temp( pid )
-                          << std::endl;
             } );
     }
 };
@@ -274,7 +268,7 @@ template <class BoundaryType, class BCTag, class ExecSpace, class Particles>
 auto createBoundaryCondition( BCTag, ExecSpace exec_space, Particles particles,
                               std::vector<BoundaryType> planes,
                               const double value,
-                              const double initial_guess = 0.5 )
+                              const double initial_guess = 1.1 )
 {
     using memory_space = typename Particles::memory_space;
     using bc_index_type = BoundaryIndexSpace<memory_space, BoundaryType>;
diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp
index 34c8cb07..afbc4f2e 100644
--- a/src/CabanaPD_Particles.hpp
+++ b/src/CabanaPD_Particles.hpp
@@ -427,13 +427,15 @@ class Particles<MemorySpace, PMB, Dimension>
         Cabana::Experimental::HDF5ParticleOutput::writeTimeStep(
             h5_config, "particles", MPI_COMM_WORLD, output_step, output_time,
             n_local, getPosition( use_reference ), sliceStrainEnergy(),
-            sliceForce(), sliceDisplacement(), sliceVelocity(), sliceDamage() );
+            sliceForce(), sliceDisplacement(), sliceVelocity(), sliceDamage(),
+            sliceTemperature() );
 #else
 #ifdef Cabana_ENABLE_SILO
         Cajita::Experimental::SiloParticleOutput::writePartialRangeTimeStep(
             "particles", local_grid->globalGrid(), output_step, output_time, 0,
             n_local, getPosition( use_reference ), sliceStrainEnergy(),
-            sliceForce(), sliceDisplacement(), sliceVelocity(), sliceDamage() );
+            sliceForce(), sliceDisplacement(), sliceVelocity(), sliceDamage(),
+            sliceTemperature() );
 #else
         log( std::cout, "No particle output enabled." );
 #endif
diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp
index b67fe7bc..2a584fbc 100644
--- a/src/CabanaPD_Solver.hpp
+++ b/src/CabanaPD_Solver.hpp
@@ -192,7 +192,7 @@ class SolverElastic
         computeEnergy( *force, *particles, *neighbors, neigh_iter_tag() );
 
         // Add boundary condition.
-        boundary_condition.apply( exec_space(), *particles );
+        boundary_condition.apply( exec_space(), *particles, 0.0 );
 
         particles->output( 0, 0.0, output_reference );
         init_time += init_timer.seconds();
@@ -231,7 +231,8 @@ class SolverElastic
             force_time += force_timer.seconds();
 
             // Add boundary condition.
-            boundary_condition.apply( exec_space(), *particles );
+            boundary_condition.apply( exec_space(), *particles,
+                                      step * inputs->timestep );
 
             // Integrate - velocity Verlet second half.
             integrate_timer.reset();