Merge pull request #572 from feathern/log_deriv2

Added finite-difference support to the log_deriv subroutine.
geodynamics · Aug 28, 2024 · ad02113 · ad02113
2 parents 919960a + e90d28d
commit ad02113
Show file tree

Hide file tree

Showing 2 changed files with 26 additions and 17 deletions.
diff --git a/src/Makefile.fdeps b/src/Makefile.fdeps
@@ -52,7 +52,7 @@ Main.o : Main.F90 Run_Parameters.o Benchmarking.o Fourier_Transform.o Timers.o S
 MakeDir.o : MakeDir.F90 
 Math_Constants.o : Math_Constants.F90 Ra_Precision.o 
 Math_Utility.o : Math_Utility.F90 
-PDE_Coefficients.o : PDE_Coefficients.F90 General_MPI.o Math_Utility.o Math_Constants.o Controls.o ProblemSize.o 
+PDE_Coefficients.o : PDE_Coefficients.F90 Finite_Difference.o General_MPI.o Math_Utility.o Math_Constants.o Controls.o ProblemSize.o 
 Parallel_Framework.o : Parallel_Framework.F90 BufferedOutput.o Structures.o Load_Balance.o General_MPI.o MPI_LAYER.o 
 Parallel_IO.o : Parallel_IO.F90 General_MPI.o Spherical_Buffer.o Legendre_Transforms.o Fourier_Transform.o MPI_LAYER.o ISendReceive.o Structures.o Parallel_Framework.o Ra_MPI_Base.o 
 ProblemSize.o : ProblemSize.F90 Ra_MPI_Base.o Timers.o BufferedOutput.o Math_Constants.o Finite_Difference.o Chebyshev_Polynomials.o Controls.o Spectral_Derivatives.o Legendre_Polynomials.o Parallel_Framework.o 

diff --git a/src/Physics/PDE_Coefficients.F90 b/src/Physics/PDE_Coefficients.F90
@@ -37,6 +37,7 @@ Module PDE_Coefficients
     Use Math_Constants
     Use Math_Utility
     Use General_MPI, Only : BCAST2D
+    Use Finite_Difference, Only : d_by_dx
     Implicit None
 
     !///////////////////////////////////////////////////////////
@@ -1583,32 +1584,40 @@ Subroutine Log_Deriv(arr1,arr2, no_log)
         ! Arr1 is assumed to be in physical space.
         ! Set no_log = .true. to take normal derivative.
 
-        Allocate(dtemp(1:n_r,1,1,2))
-        Allocate(dtemp2(1:n_r,1,1,2))
+        If (chebyshev) Then
+
+            Allocate(dtemp(1:n_r,1,1,2))
+            Allocate(dtemp2(1:n_r,1,1,2))
 
-        dtemp(:,:,:,:) = 0.0d0
-        dtemp2(:,:,:,:) = 0.0d0
-        dtemp(1:n_r,1,1,1) = arr1(1:n_r)
+            dtemp(:,:,:,:) = 0.0d0
+            dtemp2(:,:,:,:) = 0.0d0
+            dtemp(1:n_r,1,1,1) = arr1(1:n_r)
 
-        ! Transform to spectral space & de-alias
-        Call gridcp%to_Spectral(dtemp,dtemp2)
-        dtemp2((n_r*2)/3:n_r,1,1,1) = 0.0d0
+            ! Transform to spectral space & de-alias
+            Call gridcp%to_Spectral(dtemp,dtemp2)
+            dtemp2((n_r*2)/3:n_r,1,1,1) = 0.0d0
 
-        ! Take the derivative & de-alias
-        Call gridcp%d_by_dr_cp(1,2,dtemp2,1)
-        dtemp2((n_r*2)/3:n_r,1,1,2) = 0.0d0 
+            ! Take the derivative & de-alias
+            Call gridcp%d_by_dr_cp(1,2,dtemp2,1)
+            dtemp2((n_r*2)/3:n_r,1,1,2) = 0.0d0 
 
-        !Transform back to physical space.
-        Call gridcp%From_Spectral(dtemp2,dtemp)
-        arr2(:) = dtemp(:,1,1,2)
+            !Transform back to physical space.
+            Call gridcp%From_Spectral(dtemp2,dtemp)
+            arr2(:) = dtemp(:,1,1,2)
+
+            DeAllocate(dtemp,dtemp2)
+
+        Else
+
+            Call d_by_dx(arr1,arr2,1)
+
+        Endif
 
         ! If desired, convert to logarithmic derivative (default)
         If (.not. present(no_log)) Then
             arr2(:) = arr2(:)/arr1(:)
         Endif
 
-        DeAllocate(dtemp,dtemp2)
-
     End Subroutine log_deriv
 
     Subroutine Restore_Reference_Defaults