Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
trixi-framework
GitHub Repository: trixi-framework/Trixi.jl
 Path: blob/main/src/callbacks_step/stepsize_dg1d.jl
5586 views
1
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).

2
# Since these FMAs can increase the performance of many numerical algorithms,

3
# we need to opt-in explicitly.

4
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.

5
@muladd begin

6
#! format: noindent

7


8
function max_dt(u, t, mesh::TreeMesh{1},

9
                constant_speed::False, equations,

10
                dg::DG, cache)

11
    # Avoid division by zero if the speed vanishes everywhere

12
    max_scaled_speed = nextfloat(zero(t))

13


14
    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)

15
        max_lambda1 = zero(max_scaled_speed)

16
        for i in eachnode(dg)

17
            u_node = get_node_vars(u, equations, dg, i, element)

18
            lambda1, = max_abs_speeds(u_node, equations)

19
            max_lambda1 = Base.max(max_lambda1, lambda1)

20
        end

21
        inv_jacobian = cache.elements.inverse_jacobian[element] # 2 / Δx

22
        # Use `Base.max` to prevent silent failures, as `max` from `@fastmath` doesn't propagate

23
        # `NaN`s properly. See https://github.com/trixi-framework/Trixi.jl/pull/2445#discussion_r2336812323

24
        max_scaled_speed = Base.max(max_scaled_speed, inv_jacobian * max_lambda1)

25
    end

26


27
    # Factor 2 cancels with 2 from `inv_jacobian`, resulting in Δx

28
    return 2 / (nnodes(dg) * max_scaled_speed)

29
end

30


31
function max_dt(u, t, mesh::TreeMesh{1},

32
                constant_diffusivity::False, equations,

33
                equations_parabolic::AbstractEquationsParabolic,

34
                dg::DG, cache)

35
    # Avoid division by zero if the diffusivity vanishes everywhere

36
    max_scaled_diffusivity = nextfloat(zero(t))

37


38
    @batch reduction=(max, max_scaled_diffusivity) for element in eachelement(dg, cache)

39
        max_diffusivity_ = zero(max_scaled_diffusivity)

40
        for i in eachnode(dg)

41
            u_node = get_node_vars(u, equations, dg, i, element)

42
            diffusivity = max_diffusivity(u_node, equations_parabolic)

43
            max_diffusivity_ = Base.max(max_diffusivity_, diffusivity)

44
        end

45
        inv_jacobian = cache.elements.inverse_jacobian[element] # 2 / Δx

46
        # Use `Base.max` to prevent silent failures, as `max` from `@fastmath` doesn't propagate

47
        # `NaN`s properly. See https://github.com/trixi-framework/Trixi.jl/pull/2445#discussion_r2336812323

48
        max_scaled_diffusivity = Base.max(max_scaled_diffusivity,

49
                                          inv_jacobian^2 * max_diffusivity_)

50
    end

51


52
    # Factor 4 cancels with 2^2 from `inv_jacobian^2`, resulting in Δx^2

53
    return 4 / (nnodes(dg) * max_scaled_diffusivity)

54
end

55


56
function max_dt(u, t, mesh::TreeMesh{1},

57
                constant_speed::True, equations,

58
                dg::DG, cache)

59
    # Avoid division by zero if the speed vanishes everywhere,

60
    # e.g. for steady-state linear advection

61
    max_scaled_speed = nextfloat(zero(t))

62


63
    max_lambda1, = max_abs_speeds(equations)

64


65
    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)

66
        inv_jacobian = cache.elements.inverse_jacobian[element] # 2 / Δx

67
        # Use `Base.max` to prevent silent failures, as `max` from `@fastmath` doesn't propagate

68
        # `NaN`s properly. See https://github.com/trixi-framework/Trixi.jl/pull/2445#discussion_r2336812323

69
        max_scaled_speed = Base.max(max_scaled_speed, inv_jacobian * max_lambda1)

70
    end

71


72
    # Factor 2 cancels with 2 from `inv_jacobian`, resulting in Δx

73
    return 2 / (nnodes(dg) * max_scaled_speed)

74
end

75


76
function max_dt(u, t, mesh::TreeMesh, # for all dimensions

77
                constant_diffusivity::True, equations,

78
                equations_parabolic::AbstractEquationsParabolic,

79
                dg::DG, cache)

80
    # Avoid division by zero if the diffusivity vanishes everywhere

81
    max_scaled_diffusivity = nextfloat(zero(t))

82


83
    diffusivity = max_diffusivity(equations_parabolic)

84


85
    @batch reduction=(max, max_scaled_diffusivity) for element in eachelement(dg, cache)

86
        inv_jacobian = cache.elements.inverse_jacobian[element] # 2 / Δx

87
        # Note: For the currently supported parabolic equations

88
        # Diffusion & Navier-Stokes, we only have one isotropic diffusivity,

89
        # so this is valid for 1D, 2D and 3D.

90
        max_scaled_diffusivity = Base.max(max_scaled_diffusivity,

91
                                          inv_jacobian^2 * diffusivity)

92
    end

93


94
    # Factor 4 cancels with 2^2 from `inv_jacobian^2`, resulting in Δx^2

95
    return 4 / (nnodes(dg) * max_scaled_diffusivity)

96
end

97


98
function max_dt(u, t, mesh::StructuredMesh{1},

99
                constant_speed::False, equations,

100
                dg::DG, cache)

101
    # Avoid division by zero if the speed vanishes everywhere

102
    max_scaled_speed = nextfloat(zero(t))

103


104
    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)

105
        max_lambda1 = zero(max_scaled_speed)

106


107
        for i in eachnode(dg)

108
            u_node = get_node_vars(u, equations, dg, i, element)

109
            lambda1, = max_abs_speeds(u_node, equations)

110


111
            inv_jacobian = cache.elements.inverse_jacobian[i, element] # 2 / Δx

112


113
            max_lambda1 = Base.max(max_lambda1, inv_jacobian * lambda1)

114
        end

115


116
        # Use `Base.max` to prevent silent failures, as `max` from `@fastmath` doesn't propagate

117
        # `NaN`s properly. See https://github.com/trixi-framework/Trixi.jl/pull/2445#discussion_r2336812323

118
        max_scaled_speed = Base.max(max_scaled_speed, max_lambda1)

119
    end

120


121
    # Factor 2 cancels with 2 from `inv_jacobian`, resulting in Δx

122
    return 2 / (nnodes(dg) * max_scaled_speed)

123
end

124


125
function max_dt(u, t, mesh::StructuredMesh{1},

126
                constant_speed::True, equations,

127
                dg::DG, cache)

128
    # Avoid division by zero if the speed vanishes everywhere,

129
    # e.g. for steady-state linear advection

130
    max_scaled_speed = nextfloat(zero(t))

131


132
    max_lambda1, = max_abs_speeds(equations)

133


134
    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)

135
        for i in eachnode(dg)

136
            inv_jacobian = cache.elements.inverse_jacobian[i, element] # 2 / Δx

137
            # Use `Base.max` to prevent silent failures, as `max` from `@fastmath` doesn't propagate

138
            # `NaN`s properly. See https://github.com/trixi-framework/Trixi.jl/pull/2445#discussion_r2336812323

139
            max_scaled_speed = Base.max(max_scaled_speed, inv_jacobian * max_lambda1)

140
        end

141
    end

142


143
    # Factor 2 cancels with 2 from `inv_jacobian`, resulting in Δx

144
    return 2 / (nnodes(dg) * max_scaled_speed)

145
end

146


147
# Note: `max_dt` is not implemented for `StructuredMesh{1}` and `equations_parabolic` since 

148
# for the `StructuredMesh` type there is no support of parabolic terms (yet), see the overview in the docs:

149
# https://trixi-framework.github.io/Trixi.jl/stable/overview/#overview-semidiscretizations

150


151
end # @muladd

152


153