Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rework rheology database #144

Merged
merged 23 commits into from
Dec 28, 2023

Conversation

albert-de-montserrat
Copy link
Member

@albert-de-montserrat albert-de-montserrat commented Dec 27, 2023

Fixes and addresses #129. The rheology database has been rework; to make it fully type stable every rheology is contained into his own function. The rheology are organised in the following not exported submodules: GeoParams.Dislocation, GeoParams.Diffusion, GeoParams.GBS, GeoParams.NonLinearPeierls, and GeoParams.Peierls. The list of available rheologies can be accessed either by tab completion from the REPL or with:

diffusion_law_list()
dislocation_law_list()
grainboundarysliding_law_list()
nonlinearpeierls_law_list()
peierls_law_list() 

The following example does not allocate anymore

 function main(Tv, ε̇ii, η)
        # Unit system
        CharDim = SI_units(length=1000m, temperature=1000C, stress=1e7Pa, viscosity=1e20Pas)
        # Configure viscosity model
        flow_nd = SetDislocationCreep(Dislocation.diabase_Caristan_1982, CharDim)
        
        # Setup up viscosity model
        a = @allocated begin
            for i in eachindex(ε̇ii)
                η[i] = compute_viscosity_εII(flow_nd, ε̇ii[i], (;T=Tv[i]))
            end
        end

        return a
    end
   
# Numerical parameters
Ncy = 100
# Allocate arrays
ε0  = 1e-4
Tv  = rand(Ncy+1)
ε̇ii = fill(ε0, Ncy+1) 
η   = zeros(Ncy+1)

@test main(Tv, ε̇ii, η) == 0 # passes

@boriskaus
Copy link
Member

This is great, thanks a lot!
Some of the tests still indicate allocations though:

Allocations: Test Failed at /home/runner/work/GeoParams.jl/GeoParams.jl/test/test_Allocations.jl:43
  Expression: main1(Tv, ε̇ii, η) == 0
   Evaluated: 11782810 == 0

@albert-de-montserrat
Copy link
Member Author

albert-de-montserrat commented Dec 28, 2023

Looks like it still allocates in 1.6 and 1.7, I don't know if we should still support those versions though...

@albert-de-montserrat
Copy link
Member Author

Problem is that in those versions the power of a unit is always unstable...

julia> @code_warntype MPa ^2
MethodInstance for ^(::Unitful.FreeUnits{(MPa,), 𝐌 𝐋⁻¹·⁰ 𝐓⁻²·⁰, nothing}, ::Int64)
  from ^(x::Unitful.FreeUnits{N, D, nothing}, y::Integer) where {N, D} in Unitful at /Users/albert/.julia/packages/Unitful/R4J37/src/units.jl:115
Static Parameters
  N = (MPa,)
  D = 𝐌 𝐋⁻¹·⁰ 𝐓⁻²·⁰
Arguments
  #self#::Core.Const(^)
  x::Core.Const(MPa)
  y::Int64
Locals
  #40::Unitful.var"#40#41"{Int64}
Body::Any
1 ─ %1  = Unitful.:(var"#40#41")::Core.Const(Unitful.var"#40#41")
│   %2  = Core.typeof(y)::Core.Const(Int64)
│   %3  = Core.apply_type(%1, %2)::Core.Const(Unitful.var"#40#41"{Int64})
│         (#40 = %new(%3, y))
│   %5  = #40::Unitful.var"#40#41"{Int64}
│   %6  = Unitful.map(%5, $(Expr(:static_parameter, 1)))::Core.PartialStruct(Tuple{Unitful.Unit{:Pascal, 𝐌 𝐋⁻¹·⁰ 𝐓⁻²·⁰}}, Any[Core.PartialStruct(Unitful.Unit{:Pascal, 𝐌 𝐋⁻¹·⁰ 𝐓⁻²·⁰}, Any[Core.Const(6), Rational{Int64}])])
│   %7  = ::Core.Const()
│   %8  = Core.apply_type(Unitful.FreeUnits, %6, %7)::Type{Unitful.FreeUnits{_A, }} where _A
│   %9  = (%8)()::Unitful.FreeUnits{_A, , nothing} where _A
│   %10 = (*)(%9)::Any
└──       return %10

@boriskaus
Copy link
Member

I’m fine with only supporting newer versions of julia, provided this is made clear in the [compat] part. Note that CI for newer versions also fails because of an unrelated issue

Copy link
Member Author

@albert-de-montserrat albert-de-montserrat left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I addressed the remaining broken tests and allocations are now conditionally tested on >1.7, so we can keep testing for 1.6 and 1.7.

@albert-de-montserrat
Copy link
Member Author

I had to drop 1.6. 1.7 is fine -which was missing in the CI-, and updated Julia compats. Now everything should be working fine.

@albert-de-montserrat
Copy link
Member Author

Closes #129

Copy link
Member

@boriskaus boriskaus left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks a lot!

@albert-de-montserrat albert-de-montserrat merged commit 7c984f7 into main Dec 28, 2023
17 of 18 checks passed
@albert-de-montserrat albert-de-montserrat deleted the bk_type_instability_dislocation_creep branch December 28, 2023 11:42
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants