From 4471ac20fc72e5939807d13dcb63c8dc60d182d8 Mon Sep 17 00:00:00 2001 From: Loren Date: Thu, 30 Nov 2023 00:22:58 -0800 Subject: [PATCH] guarded against dividing by N^2=0 --- src/Physics/PDE_Coefficients.F90 | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/Physics/PDE_Coefficients.F90 b/src/Physics/PDE_Coefficients.F90 index 34da3d7f..7586ecf1 100644 --- a/src/Physics/PDE_Coefficients.F90 +++ b/src/Physics/PDE_Coefficients.F90 @@ -720,7 +720,9 @@ Subroutine Polytropic_ReferenceND_General() ! That's fine up to a multiplicative constant which we determine below nsquared = gravity*ref%dsdr ! N^2 (non-dimensional) up to multiplicative constant - nsquared = nsquared/nsquared(N_R) ! N^2 (non-dimensional) now correct if ND_Inner_Radius .eq. .True. + If (.not. adiabatic_polytrope) Then ! don't divide by zero! + nsquared = nsquared/nsquared(N_R) ! N^2 (non-dimensional) now correct if ND_Inner_Radius .eq. .True. + Endif ! calculate "basal" dissipation number explicitly Dissipation_Number = (poly_n + 1.0d0)/(poly_n_ad + 1.0d0) * (1.0d0/aspect_ratio) * & @@ -743,7 +745,9 @@ Subroutine Polytropic_ReferenceND_General() ref%density = ref%density/ref%density(1) ref%temperature = ref%temperature/ref%temperature(1) gravity = gravity/gravity(1) - nsquared = nsquared/nsquared(1) + If (.not. adiabatic_polytrope) Then ! don't divide by zero! + nsquared = nsquared/nsquared(1) + Endif ElseIf (ND_Volume_Average) Then ! This is the default Call Integrate_in_radius(ref%density,norm) norm = four_pi*norm/shell_volume @@ -757,13 +761,11 @@ Subroutine Polytropic_ReferenceND_General() norm = four_pi*norm/shell_volume gravity = gravity/norm - Call Integrate_in_radius(nsquared,norm) - norm = four_pi*norm/shell_volume - nsquared = nsquared/norm - - Call Integrate_in_radius(nsquared,norm) - norm = four_pi*norm/shell_volume - nsquared = nsquared/norm + If (.not. adiabatic_polytrope) Then ! don't divide by zero! + Call Integrate_in_radius(nsquared,norm) + norm = four_pi*norm/shell_volume + nsquared = nsquared/norm + Endif ! ref%temperature(N_R) is now T_inner/volav(T); gravity(N_R) = g_inner/volav(g) Dissipation_Number = Dissipation_Number*ref%temperature(N_R)/gravity(N_R)