revision of comments

PTsolvers · Dec 20, 2023 · 6f2f524 · 6f2f524
1 parent 29f5692
commit 6f2f524
Show file tree

Hide file tree

Showing 3 changed files with 11 additions and 17 deletions.
diff --git a/miniapps/benchmarks/stokes2D/ShearHeating/Shearheating2D.jl b/miniapps/benchmarks/stokes2D/ShearHeating/Shearheating2D.jl
@@ -1,8 +1,5 @@
 # Benchmark of Duretz et al. 2014
 # http://dx.doi.org/10.1002/2014GL060438
-using CUDA
-CUDA.allowscalar(false) # for safety
-
 using JustRelax, JustRelax.DataIO, JustPIC
 import JustRelax.@cell
 
@@ -153,10 +150,10 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", save_vtk =false)
 
     # Buoyancy forces
     ρg               = @zeros(ni...), @zeros(ni...)
-    for _ in 1:1
-        @parallel (JustRelax.@idx ni) compute_ρg!(ρg[2], phase_ratios.center, rheology, (T=thermal.Tc, P=stokes.P))
-        @parallel init_P!(stokes.P, ρg[2], xci[2])
-    end
+
+    @parallel (JustRelax.@idx ni) compute_ρg!(ρg[2], phase_ratios.center, rheology, (T=thermal.Tc, P=stokes.P))
+    @parallel init_P!(stokes.P, ρg[2], xci[2])
+
     # Rheology
     η                = @zeros(ni...)
     args             = (; T = thermal.Tc, P = stokes.P, dt = Inf)
@@ -177,7 +174,6 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", save_vtk =false)
     )
     ## Compression and not extension - fix this
     εbg              = 5e-14
-    # εbg              = 1e-16
     stokes.V.Vx .= PTArray([ -(x - lx/2) * εbg for x in xvi[1], _ in 1:ny+2])
     stokes.V.Vy .= PTArray([ (ly - abs(y)) * εbg for _ in 1:nx+2, y in xvi[2]])
     flow_bcs!(stokes, flow_bcs) # apply boundary conditions
@@ -258,7 +254,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", save_vtk =false)
           @views thermal.T[2:end-1, :] .= T_buffer
           temperature2center!(thermal)
 
-          @parallel (@idx ni) compute_SH!(
+          @parallel (@idx ni) compute_shear_heating(
             thermal.shear_heating,
             @tensor_center(stokes.τ),
             @tensor_center(stokes.τ_o),

diff --git a/miniapps/benchmarks/stokes2D/ShearHeating/Shearheating_rheology.jl b/miniapps/benchmarks/stokes2D/ShearHeating/Shearheating_rheology.jl
@@ -2,9 +2,7 @@
 
 function init_rheologies(; is_TP_Conductivity=true)
 
-    # Dislocation and Diffusion creep
-    # Matrix            = LinearViscous(η=1e23)
-    # Inclusion         = LinearViscous(η=1e22)
+    # Dislocation creep
     Matrix        = DislocationCreep(A=3.20e-20, n=3.0, E=276e3, V=0e0,  r=0.0, R=8.3145)
     Inclusion     = DislocationCreep(A=3.16e-26, n=3.3, E=186e3, V=0e0,  r=0.0, R=8.3145)
 
@@ -17,7 +15,7 @@ function init_rheologies(; is_TP_Conductivity=true)
             d = 0.0,
         )
     else
-        K_Matrix = ConstantConductivity(k = 2.5)
+        ConstantConductivity(k = 2.5)
     end
 
     K_Inclusion = if is_TP_Conductivity

diff --git a/src/thermal_diffusion/Shearheating.jl b/src/thermal_diffusion/Shearheating.jl
@@ -1,14 +1,14 @@
-@parallel_indices (I...) function compute_SH!(
+@parallel_indices (I...) function compute_shear_heating(
     shear_heating, τ::NTuple{N,T}, τ_old::NTuple{N,T}, ε::NTuple{N,T}, rheology, dt
 ) where {N,T}
     _Gdt = inv(fn_ratio(get_G, rheology) * dt)
-    τij, τij_o, εij = JustRelax.cache_tensors(τ, τ_old, ε, I...)
+    τij, τij_o, εij = cache_tensors(τ, τ_old, ε, I...)
     εij_el = @. 0.5 * ((τij - τij_o) * _Gdt)
     shear_heating[I...] = compute_shearheating(rheology, τij, εij, εij_el)
     return nothing
 end
 
-@parallel_indices (I...) function compute_SH!(
+@parallel_indices (I...) function compute_shear_heating(
     shear_heating,
     τ::NTuple{N,T},
     τ_old::NTuple{N,T},
@@ -19,7 +19,7 @@ end
 ) where {N,T}
     phase = @inbounds phase_ratios[I...]
     _Gdt = inv(fn_ratio(get_G, rheology, phase) * dt)
-    τij, τij_o, εij = JustRelax.cache_tensors(τ, τ_old, ε, I...)
+    τij, τij_o, εij = cache_tensors(τ, τ_old, ε, I...)
     εij_el = @. 0.5 * ((τij - τij_o) * _Gdt)
     shear_heating[I...] = compute_shearheating(rheology, phase, τij, εij, εij_el)
     return nothing