Compute melt fraction (#223)

* add `compute_melt_fraction` * fix melting.jl and miniapps; add test * fix test; format * fix tests; address requested changes
PTsolvers · Aug 26, 2024 · f91fd23 · f91fd23
1 parent 3d5df73
commit f91fd23
Show file tree

Hide file tree

Showing 8 changed files with 538 additions and 3 deletions.
diff --git a/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim.jl b/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim.jl
@@ -147,8 +147,6 @@ function init_rheology(CharDim; is_compressible = false, steady_state=true)
         creep_rock  = DislocationCreep(; A=1.67e-24, n=3.5, E=1.87e5, V=6e-6, r=0.0, R=8.3145)
         creep_magma = DislocationCreep(; A=1.67e-24, n=3.5, E=1.87e5, V=6e-6, r=0.0, R=8.3145)
         creep_air   = LinearViscous(; η=1e18 * Pa * s)
-        β_rock      = 6.0e-11
-        β_magma     = 6.0e-11
     end
     g = 9.81m/s^2
     rheology = (
@@ -453,8 +451,8 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false)
             checkpointing_hdf5(figdir, stokes, thermal.T, t, dt)
 
             if igg.me == 0
+                velocity2vertex!(Vx_v, Vy_v, @velocity(stokes)...)
                 if do_vtk
-                    velocity2vertex!(Vx_v, Vy_v, @velocity(stokes)...)
                     data_v = (;
                         T   = Array(ustrip.(dimensionalize(thermal.T[2:(end - 1), :], C, CharDim))),
                         τxy = Array(ustrip.(dimensionalize(stokes.τ.xy, s^-1, CharDim))),

diff --git a/src/common.jl b/src/common.jl
@@ -54,6 +54,9 @@ export compute_ρg!
 include("rheology/Viscosity.jl")
 export compute_viscosity!
 
+include("rheology/Melting.jl")
+export compute_melt_fraction!
+
 # include("thermal_diffusion/DiffusionExplicit.jl")
 # export ThermalParameters
 

diff --git a/src/ext/AMDGPU/2D.jl b/src/ext/AMDGPU/2D.jl
@@ -226,6 +226,17 @@ function JR2D.compute_ρg!(ρg::ROCArray, phase_ratios::JustRelax.PhaseRatio, rh
     return compute_ρg!(ρg, phase_ratios, rheology, args)
 end
 
+## Melt fraction
+function JR2D.compute_melt_fraction!(ϕ::ROCArray, rheology, args)
+    return compute_melt_fraction!(ϕ, rheology, args)
+end
+
+function JR2D.compute_melt_fraction!(
+    ϕ::ROCArray, phase_ratios::JustRelax.PhaseRatio, rheology, args
+)
+    return compute_melt_fraction!(ϕ, phase_ratios, rheology, args)
+end
+
 # Interpolations
 function JR2D.temperature2center!(::AMDGPUBackendTrait, thermal::JustRelax.ThermalArrays)
     return _temperature2center!(thermal)

diff --git a/src/ext/AMDGPU/3D.jl b/src/ext/AMDGPU/3D.jl
@@ -226,6 +226,17 @@ function JR3D.compute_ρg!(ρg::ROCArray, phase_ratios::JustRelax.PhaseRatio, rh
     return compute_ρg!(ρg, phase_ratios, rheology, args)
 end
 
+## Melt fraction
+function JR3D.compute_melt_fraction!(ϕ::ROCArray, rheology, args)
+    return compute_melt_fraction!(ϕ, rheology, args)
+end
+
+function JR3D.compute_melt_fraction!(
+    ϕ::ROCArray, phase_ratios::JustRelax.PhaseRatio, rheology, args
+)
+    return compute_melt_fraction!(ϕ, phase_ratios, rheology, args)
+end
+
 # Interpolations
 function JR3D.temperature2center!(::AMDGPUBackendTrait, thermal::JustRelax.ThermalArrays)
     return _temperature2center!(thermal)

diff --git a/src/ext/CUDA/2D.jl b/src/ext/CUDA/2D.jl
@@ -222,6 +222,17 @@ function JR2D.compute_ρg!(ρg::CuArray, phase_ratios::JustRelax.PhaseRatio, rhe
     return compute_ρg!(ρg, phase_ratios, rheology, args)
 end
 
+## Melt fraction
+function JR2D.compute_melt_fraction!(ϕ::CuArray, rheology, args)
+    return compute_melt_fraction!(ϕ, rheology, args)
+end
+
+function JR2D.compute_melt_fraction!(
+    ϕ::CuArray, phase_ratios::JustRelax.PhaseRatio, rheology, args
+)
+    return compute_melt_fraction!(ϕ, phase_ratios, rheology, args)
+end
+
 # Interpolations
 function JR2D.temperature2center!(::CUDABackendTrait, thermal::JustRelax.ThermalArrays)
     return _temperature2center!(thermal)

diff --git a/src/ext/CUDA/3D.jl b/src/ext/CUDA/3D.jl
@@ -233,6 +233,17 @@ function JR3D.compute_ρg!(ρg::CuArray, phase_ratios::JustRelax.PhaseRatio, rhe
     return compute_ρg!(ρg, phase_ratios, rheology, args)
 end
 
+## Melt fraction
+function JR3D.compute_melt_fraction!(ϕ::CuArray, rheology, args)
+    return compute_melt_fraction!(ϕ, rheology, args)
+end
+
+function JR3D.compute_melt_fraction!(
+    ϕ::CuArray, phase_ratios::JustRelax.PhaseRatio, rheology, args
+)
+    return compute_melt_fraction!(ϕ, phase_ratios, rheology, args)
+end
+
 # Interpolations
 function JR3D.temperature2center!(::CUDABackendTrait, thermal::JustRelax.ThermalArrays)
     return _temperature2center!(thermal)

diff --git a/src/rheology/Melting.jl b/src/rheology/Melting.jl
@@ -0,0 +1,31 @@
+function compute_melt_fraction!(ϕ, rheology, args)
+    ni = size(ϕ)
+    @parallel (@idx ni) compute_melt_fraction_kernel!(ϕ, rheology, args)
+end
+
+@parallel_indices (I...) function compute_melt_fraction_kernel!(ϕ, rheology, args)
+    args_ijk = ntuple_idx(args, I...)
+    ϕ[I...] = compute_melt_frac(rheology, args_ijk)
+    return nothing
+end
+
+@inline function compute_melt_frac(rheology, args)
+    return compute_meltfraction(rheology, args)
+end
+
+function compute_melt_fraction!(ϕ, phase_ratios, rheology, args)
+    ni = size(ϕ)
+    @parallel (@idx ni) compute_melt_fraction_kernel!(ϕ, phase_ratios, rheology, args)
+end
+
+@parallel_indices (I...) function compute_melt_fraction_kernel!(
+    ϕ, phase_ratios, rheology, args
+)
+    args_ijk = ntuple_idx(args, I...)
+    ϕ[I...] = compute_melt_frac(rheology, args_ijk, phase_ratios[I...])
+    return nothing
+end
+
+@inline function compute_melt_frac(rheology, args, phase_ratios)
+    return compute_meltfraction_ratio(phase_ratios, rheology, args)
+end