add visc and buoyancy calculations inside the solver

PTsolvers · May 3, 2024 · ae5b77f · ae5b77f
1 parent b286112
commit ae5b77f
Show file tree

Hide file tree

Showing 2 changed files with 16 additions and 0 deletions.
diff --git a/src/stokes/Stokes2D.jl b/src/stokes/Stokes2D.jl
@@ -320,6 +320,10 @@ function _solve!(
     θ = @zeros(ni...)
     Vx_on_Vy = @zeros(size(stokes.V.Vy))
 
+    # compute buoyancy forces and viscosity
+    compute_ρg!(ρg[end], phase_ratios, rheology, args)
+    compute_viscosity!(stokes, phase_ratios, args, rheology, cutoff_visc)
+
     while iter < 2 || (err > ϵ && iter ≤ iterMax)
         wtime0 += @elapsed begin
             @parallel (@idx ni) compute_∇V!(stokes.∇V, @velocity(stokes)..., _di...)
@@ -502,6 +506,10 @@ function _solve!(
     end
     Vx_on_Vy = @zeros(size(stokes.V.Vy))
 
+    # compute buoyancy forces and viscosity
+    compute_ρg!(ρg[end], phase_ratios, rheology, args)
+    compute_viscosity!(stokes, phase_ratios, args, rheology, cutoff_visc)
+
     while iter ≤ iterMax
         iterMin < iter && err < ϵ && break
 

diff --git a/src/stokes/Stokes3D.jl b/src/stokes/Stokes3D.jl
@@ -202,6 +202,10 @@ function _solve!(
     λ = @zeros(ni...)
     θ = @zeros(ni...)
 
+    # compute buoyancy forces and viscosity
+    compute_ρg!(ρg[end], phase_ratios, rheology, args)
+    compute_viscosity!(stokes, phase_ratios, args, rheology, cutoff_visc)
+
     # solver loop
     wtime0 = 0.0
     while iter < 2 || (err > ϵ && iter ≤ iterMax)
@@ -382,6 +386,10 @@ function _solve!(
     wtime0 = 0.0
     ητ = deepcopy(η)
 
+    # compute buoyancy forces and viscosity
+    compute_ρg!(ρg[end], phase_ratios, rheology, args)
+    compute_viscosity!(stokes, phase_ratios, args, rheology, cutoff_visc)
+
     while iter < 2 || (err > ϵ && iter ≤ iterMax)
         wtime0 += @elapsed begin
             # ~preconditioner