Clean up treatment of floating point precision (#346)

Author: glwagner

src/Bathymetry.jl

@@ -165,7 +165,7 @@ function regrid_bathymetry(target_grid;
     Nyn = length(φ_data)
     Nzn = 1
 
-    native_grid = LatitudeLongitudeGrid(arch;
+    native_grid = LatitudeLongitudeGrid(arch, Float32;
                                         size = (Nxn, Nyn, Nzn),
                                         latitude  = (φ₁_data, φ₂_data),
                                         longitude = (λ₁_data, λ₂_data),
@@ -234,7 +234,7 @@ function interpolate_bathymetry_in_passes(native_z, target_grid;
 
         @debug "Bathymetry interpolation pass $pass with size $new_size"
 
-        new_grid = LatitudeLongitudeGrid(architecture(target_grid),
+        new_grid = LatitudeLongitudeGrid(architecture(target_grid), Float32,
                                          size = new_size,
                                          latitude = (latitude[1],  latitude[2]),
                                          longitude = (longitude[1], longitude[2]),

src/DataWrangling/ECCO/ECCO.jl

@@ -115,7 +115,7 @@ function empty_ECCO_field(metadata::ECCOMetadata;
         sz = (Nx, Ny)
     end
 
-    grid = LatitudeLongitudeGrid(architecture; halo, longitude, latitude, z,
+    grid = LatitudeLongitudeGrid(architecture, Float32; halo, longitude, latitude, z,
                                  size = sz,
                                  topology = (TX, TY, TZ))
 

src/OceanSeaIceModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl

@@ -12,7 +12,7 @@ end
 convert_if_number(FT, a::Number) = convert(FT, a)
 convert_if_number(FT, a) = a
 
-function CoefficientBasedFluxes(FT = Float64;
+function CoefficientBasedFluxes(FT = Oceananigans.defaults.FloatType;
                                 drag_coefficient = 1e-3,
                                 gravitational_acceleration = g_Earth,
                                 heat_transfer_coefficient = drag_coefficient,

src/OceanSeaIceModels/InterfaceComputations/interface_states.jl

@@ -59,7 +59,7 @@ struct WaterMoleFraction{FT, C}
     salinity_constituents :: C
 end
 
-function WaterMoleFraction(FT=Float64)
+function WaterMoleFraction(FT=Oceananigans.defaults.FloatType)
     water_molar_mass = convert(FT, 18.02)
 
     # TODO: find reference for these

src/OceanSeaIceModels/InterfaceComputations/radiation.jl

@@ -38,7 +38,7 @@ Keyword Arguments
 - `sea_ice_albedo`: The albedo of the sea ice surface. Default: `0.7`.
 - `stefan_boltzmann_constant`: The Stefan-Boltzmann constant. Default: `5.67e-8`.
 """
-function Radiation(arch = CPU(), FT=Float64;
+function Radiation(arch = CPU(), FT=Oceananigans.defaults.FloatType;
                    ocean_emissivity = 0.97,
                    sea_ice_emissivity = 1.0,
                    ocean_albedo = LatitudeDependentAlbedo(FT),

src/OceanSeaIceModels/InterfaceComputations/roughness_lengths.jl

@@ -13,7 +13,7 @@ struct ScalarRoughnessLength{FT, V, R}
 end
 
 """
-    ScalarRoughnessLength([FT=Float64];
+    ScalarRoughnessLength(FT = Float64;
                           air_kinematic_viscosity = temperature_dependent_viscosity,
                           reynolds_number_scaling_function = empirical_scaling_function,
                           maximum_roughness_length = 1.6e-4)
@@ -28,7 +28,7 @@ Keyword Arguments
 - `reynolds_number_scaling_function::Function`: The function to compute the Reynolds number scaling factor.
 - `maximum_roughness_length::Float`: The maximum roughness length value. Defaults to `1.6e-4`.
 """
-function ScalarRoughnessLength(FT=Float64;
+function ScalarRoughnessLength(FT=Oceananigans.defaults.FloatType;
                                air_kinematic_viscosity = TemperatureDependentAirViscosity(FT),
                                reynolds_number_scaling_function = ReynoldsScalingFunction(FT),
                                maximum_roughness_length = 1.6e-4) # Values from COARE3.6
@@ -39,7 +39,7 @@ function ScalarRoughnessLength(FT=Float64;
 end
 
 """
-    MomentumRoughnessLength([FT=Float64];
+    MomentumRoughnessLength(FT = Float64;
                             gravitational_acceleration = default_gravitational_acceleration,
                             maximum_roughness_length = 1.0,
                             air_kinematic_viscosity = TemperatureDependentAirViscosity(FT),
@@ -58,7 +58,7 @@ Keyword Arguments
 - `gravity_wave_parameter`: The wave parameter. Default: 0.011.
 - `laminar_parameter`: The laminar parameter. Default: 0.11.
 """
-function MomentumRoughnessLength(FT=Float64;
+function MomentumRoughnessLength(FT=Oceananigans.defaults.FloatType;
                                  gravitational_acceleration = default_gravitational_acceleration,
                                  maximum_roughness_length = 1.0, # An estimate?
                                  air_kinematic_viscosity = TemperatureDependentAirViscosity(FT),
@@ -72,7 +72,7 @@ function MomentumRoughnessLength(FT=Float64;
                                           convert(FT, maximum_roughness_length))
 end
 
-function default_roughness_lengths(FT=Float64)
+function default_roughness_lengths(FT=Oceananigans.defaults.FloatType)
     momentum    = MomentumRoughnessLength(FT)
     temperature = ScalarRoughnessLength(FT)
     water_vapor = ScalarRoughnessLength(FT)
@@ -88,7 +88,7 @@ struct TemperatureDependentAirViscosity{FT}
 end
 
 """
-    TemperatureDependentAirViscosity([FT = Float64;
+    TemperatureDependentAirViscosity([FT = Oceananigans.defaults.FloatType;
                                       C₀ = 1.326e-5,
                                       C₁ = C₀ * 6.542e-3,
                                       C₂ = C₀ * 8.301e-6,
@@ -100,7 +100,7 @@ viscosity of air as
 C₀ + C₁ T + C₂ T^2 + C₃ T^3.
 ```
 """
-function TemperatureDependentAirViscosity(FT = Float64;
+function TemperatureDependentAirViscosity(FT = Oceananigans.defaults.FloatType;
                                           C₀ = 1.326e-5,
                                           C₁ = C₀ * 6.542e-3,
                                           C₂ = C₀ * 8.301e-6,
@@ -148,15 +148,15 @@ struct ReynoldsScalingFunction{FT}
 end
 
 """
-    ReynoldsScalingFunction(FT = Float64; A = 5.85e-5, b = 0.72)
+    ReynoldsScalingFunction(FT=Float64; A=5.85e-5, b=0.72)
 
 Empirical fit of the scalar roughness length with roughness Reynolds number `R★ = u★ ℓu / ν`.
 Edson et al. (2013), equation (28).
 ```math
     ℓs = A / R★ ^ b
 ```
 """
-ReynoldsScalingFunction(FT = Float64; A = 5.85e-5, b = 0.72) = 
+ReynoldsScalingFunction(FT = Oceananigans.defaults.FloatType; A = 5.85e-5, b = 0.72) = 
     ReynoldsScalingFunction(convert(FT, A), convert(FT, b))
 
 @inline (s::ReynoldsScalingFunction)(R★, args...) = ifelse(R★ == 0, convert(eltype(R★), 0), s.A / R★ ^ s.b)

src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl

@@ -106,7 +106,7 @@ Keyword Arguments
 - `solver_tolerance`: The tolerance for convergence. Default: 1e-8.
 - `solver_maxiter`: The maximum number of iterations. Default: 100.
 """
-function SimilarityTheoryFluxes(FT::DataType = Float64;
+function SimilarityTheoryFluxes(FT::DataType = Oceananigans.defaults.FloatType;
                                 gravitational_acceleration = g_Earth,
                                 von_karman_constant = 0.4,
                                 turbulent_prandtl_number = 1,
@@ -462,7 +462,7 @@ end
 end
 
 # Edson et al. (2013)
-function edson_stability_functions(FT = Float64)
+function edson_stability_functions(FT=Oceananigans.defaults.FloatType)
     ψu = EdsonMomentumStabilityFunction{FT}()
     ψc = EdsonScalarStabilityFunction{FT}()
     return SimilarityScales(ψu, ψc, ψc)
@@ -559,7 +559,7 @@ end
     return ifelse(stable, Ψ_stable, Ψ_unstable)
 end
 
-function atmosphere_sea_ice_stability_functions(FT=Float64)
+function atmosphere_sea_ice_stability_functions(FT=Oceananigans.defaults.FloatType)
     stable_momentum = PaulsonMomentumStabilityFunction{FT}()
     unstable_momentum = ShebaMomentumStabilityFunction{FT}()
     momentum = SplitStabilityFunction(stable_momentum, unstable_momentum)

src/OceanSeaIceModels/InterfaceComputations/tabulated_albedo.jl

@@ -79,7 +79,7 @@ Keyword Arguments
 - `φ_values`: The latitude values for the table. Default: `(0:2:90) ./ 180 * π`.
 - `𝓉_values`: The transmissivity values for the table. Default: `0:0.05:1`.
 """
-function TabulatedAlbedo(arch = CPU(), FT = Float64;
+function TabulatedAlbedo(arch = CPU(), FT = Oceananigans.defaults.FloatType;
                          S₀ = convert(FT, 1365),
                          α_table  = α_payne,
                          φ_values = (0:2:90) ./ 180 * π,

src/OceanSeaIceModels/PrescribedAtmospheres.jl

@@ -63,9 +63,10 @@ end
 Base.show(io::IO, p::ConstitutiveParameters) = print(io, summary(p))
 
 """
-    ConstitutiveParameters(FT; gas_constant       = 8.3144598,
-                               dry_air_molar_mass = 0.02897,
-                               water_molar_mass   = 0.018015)
+    ConstitutiveParameters(FT = Float64;
+                           gas_constant       = 8.3144598,
+                           dry_air_molar_mass = 0.02897,
+                           water_molar_mass   = 0.018015)
 
 Construct a set of parameters that define the density of moist air,
 
@@ -85,7 +86,7 @@ where
 
 For more information see [reference docs].
 """
-function ConstitutiveParameters(FT = Float64;
+function ConstitutiveParameters(FT = Oceananigans.defaults.FloatType;
                                 gas_constant       = 8.3144598,
                                 dry_air_molar_mass = 0.02897,
                                 water_molar_mass   = 0.018015)
@@ -130,7 +131,7 @@ Base.show(io::IO, p::HeatCapacityParameters) = print(io, summary(p))
 
 Isobaric heat capacities.
 """
-function HeatCapacityParameters(FT = Float64;
+function HeatCapacityParameters(FT = Oceananigans.defaults.FloatType;
                                 dry_air_adiabatic_exponent = 2/7,
                                 water_vapor_heat_capacity = 1859,
                                 liquid_water_heat_capacity = 4181,
@@ -173,7 +174,7 @@ end
 
 Base.show(io::IO, p::PhaseTransitionParameters) = print(io, summary(p))
 
-function PhaseTransitionParameters(FT = Float64;
+function PhaseTransitionParameters(FT = Oceananigans.defaults.FloatType;
                                    reference_vaporization_enthalpy = 2500800,
                                    reference_sublimation_enthalpy = 2834400,
                                    reference_temperature = 273.16,
@@ -244,7 +245,7 @@ function Base.show(io::IO, p::PrescribedAtmosphereThermodynamicsParameters)
         "    └── total_ice_nucleation_temperature (Tⁱ): ", prettysummary(pt.total_ice_nucleation_temperature))
 end
 
-function PrescribedAtmosphereThermodynamicsParameters(FT=Float64;
+function PrescribedAtmosphereThermodynamicsParameters(FT = Oceananigans.defaults.FloatType;
                                                       constitutive = ConstitutiveParameters(FT),
                                                       phase_transitions = PhaseTransitionParameters(FT),
                                                       heat_capacity = HeatCapacityParameters(FT))

src/OceanSeaIceModels/freezing_limited_ocean_temperature.jl

@@ -23,7 +23,8 @@ the `T < Tₘ` except for heating to allow temperature to increase.
 
 The melting temperature is a function of salinity and is controlled by the `liquidus`.
 """
-FreezingLimitedOceanTemperature(FT::DataType=Float64) = FreezingLimitedOceanTemperature(LinearLiquidus(FT))
+FreezingLimitedOceanTemperature(FT::DataType=Oceananigans.defaults.FloatType) =
+    FreezingLimitedOceanTemperature(LinearLiquidus(FT))
 
 const FreezingLimitedCoupledModel = OceanSeaIceModel{<:FreezingLimitedOceanTemperature}