Merge pull request #566 from feathern/fd_second_deriv

Bug Fix:  Diagnostics and Second Radial Derivatives when Using Finite-Difference Mode

CHANGELOG.md

@@ -50,6 +50,8 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm
 
 - Running the configure command now requires the correct combination of library locations to be specified; \[Brandon Lazard; 6-21-2024] [#535](https://github.com/geodynamics/Rayleigh/pull/535) 
 
+- Fixed a bug that was causing Rayleigh to crash during output if finite-difference mode was active and second-order derivatives in radius were required for output. \[Nick Featherstone; 6-26-2024; [#566](https://github.com/geodynamics/Rayleigh/pull/566) \]
+
 
 ## [1.2.0] - 5-29-2024
 

src/Diagnostics/Diagnostics_Second_Derivatives.F90

@@ -27,6 +27,7 @@ Module Diagnostics_Second_Derivatives
     Use Diagnostics_Base
     Use Structures
     Use Spectral_Derivatives
+    Use Finite_Difference, Only : d_by_dx3d3
     Implicit None
 
 
@@ -378,20 +379,31 @@ Subroutine Compute_Second_Derivatives(inbuffer)
         !Move to p1b configuration
         Call d2buffer%reform()
         Call d2buffer%construct('p1a')
-        Call gridcp%To_Spectral(d2buffer%p1b,d2buffer%p1a)
-        Call gridcp%dealias_buffer(d2buffer%p1a)
+        If (chebyshev) Then
+            Call gridcp%To_Spectral(d2buffer%p1b,d2buffer%p1a)
+            Call gridcp%dealias_buffer(d2buffer%p1a)
+        Else
+            d2buffer%p1a = d2buffer%p1b
+        Endif
         d2buffer%p1b = 0.0
         d2buffer%config='p1a'
 
         !Write(6,*)'Steps 4 complete.'
 
         !////////////////////////////////////////////////////////////
         ! Steps 5-6:  Load d2_by_dr2 and d2_by_drdt into the buffer
-        Do i = 1, nddfields*2
-            j = i+nddfields*2
-            Call gridcp%d_by_dr_cp(i,j,d2buffer%p1a,1)
-        Enddo
-        Call gridcp%From_Spectral(d2buffer%p1a,d2buffer%p1b)
+        If (chebyshev) Then
+            Do i = 1, nddfields*2
+                j = i+nddfields*2
+                Call gridcp%d_by_dr_cp(i,j,d2buffer%p1a,1)
+            Enddo
+            Call gridcp%From_Spectral(d2buffer%p1a,d2buffer%p1b)
+        Else
+            Do i = 1, nddfields*2
+                j = i+nddfields*2
+                Call d_by_dx3d3(i,j,d2buffer%p1a,1)
+            Enddo
+        Endif
         d2buffer%p1a=d2buffer%p1b
         Call d2buffer%deconstruct('p1b')
 

src/Makefile.fdeps

@@ -25,7 +25,7 @@ Diagnostics_Mean_Correction.o : Diagnostics_Mean_Correction.F90 indices.F Diagno
 Diagnostics_Miscellaneous.o : Diagnostics_Miscellaneous.F90 indices.F Diagnostics_Base.o 
 Diagnostics_Poynting_Flux.o : Diagnostics_Poynting_Flux.F90 indices.F Diagnostics_Base.o 
 Diagnostics_Scalars.o : Diagnostics_Scalars.F90 indices.F Controls.o Diagnostics_Base.o 
-Diagnostics_Second_Derivatives.o : Diagnostics_Second_Derivatives.F90 indices.F Spectral_Derivatives.o Structures.o Diagnostics_Base.o 
+Diagnostics_Second_Derivatives.o : Diagnostics_Second_Derivatives.F90 indices.F Finite_Difference.o Spectral_Derivatives.o Structures.o Diagnostics_Base.o 
 Diagnostics_Thermal_Energies.o : Diagnostics_Thermal_Energies.F90 indices.F Diagnostics_ADotGradB.o Diagnostics_Base.o 
 Diagnostics_Thermal_Equation.o : Diagnostics_Thermal_Equation.F90 indices.F Diagnostics_ADotGradB.o Diagnostics_Base.o 
 Diagnostics_Thermodynamic_Gradients.o : Diagnostics_Thermodynamic_Gradients.F90 indices.F Diagnostics_Base.o