update thermo source for spherical2d geometry, this adds -p divU sour…

…ce term to internal energy
zhichen3 · Sep 23, 2024 · 8fe9d89 · 8fe9d89
1 parent fb34ee1
commit 8fe9d89
Showing 1 changed file with 12 additions and 2 deletions.
diff --git a/Source/sources/Castro_thermo.cpp b/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) -