update thermo source for spherical2d geometry, this adds -p divU source term to internal energy (#2965)

This updates fill_thermo_source for spherical2d geometry, even though this seems to only affect sdc instead of ctu.

Author: zhichen3

Diff for: Source/sources/Castro_thermo.cpp

@@ -169,8 +169,18 @@ Castro::fill_thermo_source (MultiFab& state_in, MultiFab& thermo_src)
       }
 
 #if AMREX_SPACEDIM >= 2
-      src(i,j,k,UEINT) += -0.5_rt*(U(i,j+1,k,UMY)/U(i,j+1,k,URHO) -
-                                   U(i,j-1,k,UMY)/U(i,j-1,k,URHO))/dx[1];
+      if (coord == 2) {
+          Real thetap = prob_lo[1] + (static_cast<Real>(j) + 1.5_rt)*dx[1];
+          Real thetam = prob_lo[1] + (static_cast<Real>(j) - 0.5_rt)*dx[1];
+          Real theta = 0.5_rt*(thetam + thetap);
+
+          src(i,j,k,UEINT) += -0.5_rt*(std::sin(thetap)*U(i,j+1,k,UMY)/U(i,j+1,k,URHO) -
+                                       std::sin(thetam)*U(i,j-1,k,UMY)/U(i,j-1,k,URHO)) /
+                                      (r*std::sin(theta)*dx[1]);
+      } else {
+          src(i,j,k,UEINT) += -0.5_rt*(U(i,j+1,k,UMY)/U(i,j+1,k,URHO) -
+                                       U(i,j-1,k,UMY)/U(i,j-1,k,URHO))/dx[1];
+      }
 #endif
 #if AMREX_SPACEDIM == 3
       src(i,j,k,UEINT) += -0.5_rt*(U(i,j,k+1,UMZ)/U(i,j,k+1,URHO) -