From cf9100b951c730e62d2b926aff9cedac81a294b3 Mon Sep 17 00:00:00 2001 From: Karthik Rishinarada Date: Wed, 5 Mar 2025 11:48:16 +0530 Subject: [PATCH 1/5] buckingham add --- src/systems/buckinghum.jl | 174 ++++++++++++++++++++++++++++++++++++++ test/buckinghum.jl | 21 +++++ 2 files changed, 195 insertions(+) create mode 100644 src/systems/buckinghum.jl create mode 100644 test/buckinghum.jl diff --git a/src/systems/buckinghum.jl b/src/systems/buckinghum.jl new file mode 100644 index 0000000000..e5341e1c1c --- /dev/null +++ b/src/systems/buckinghum.jl @@ -0,0 +1,174 @@ +using ModelingToolkit +using DynamicQuantities; +using Unitful +using LinearAlgebra + + +""" + no_of_fundamental_dims(mp::Vector{Any}) + Finds number of fundamental dimensions +""" +function no_of_fundamental_dims(mp) + fundamental_dimensions = 0 + for val in mp + if val == 1 + fundamental_dimensions += 1 + end + end + + return fundamental_dimensions +end + + +""" + get_dims_of_vars(Vector{Any}, Number, Vector{Any}) + + Gets units of each variable in an array. Returns a matrix where each row corresponds + to binary + Example : kg*m^3 + Each index corresponds to length, mass, time, current, luminosity and temperature + respectively. +""" +function get_dims_of_vars(dims_vars, total_vars, dim_map) + # For every single variable it contains row of 1s and 0s mentioning which unit is present + dims_of_all_vars = zeros(total_vars, 6) + + for (ind, dim) in enumerate(dims_vars) + temp_dims = 0 + temp_dims_arr = zeros(6) + if ulength(dim) != 0 + dim_map[1] = 1 + temp_dims_arr[1] = ulength(dim) + temp_dims += 1 + end + if umass(dim) != 0 + dim_map[2] = 1 + temp_dims_arr[2] = umass(dim) + temp_dims += 1 + end + if utime(dim) != 0 + dim_map[3] = 1 + temp_dims_arr[3] = utime(dim) + temp_dims +=1 + end + if ucurrent(dim) != 0 + dim_map[4] = 1 + temp_dims_arr[4] = ucurrent(dim) + temp_dims +=1 + end + if utemperature(dim) != 0 + dim_map[5] = 1 + temp_dims_arr[5] = utemperature(dim) + + temp_dims +=1 + end + if uluminosity(dim) != 0 + dim_map[6] = 1 + temp_dims_arr[6] = uluminosity(dim) + temp_dims +=1 + end + dims_of_all_vars[ind,:] = temp_dims_arr' + end + + return dims_of_all_vars +end + + +""" + # find_pi_term_exponents(Matrix{Float64}, Number, Number) + + Finds PI terms and returns the exponents and indices of repeating variables and + non repeating index in the form of a dictionary +""" +function find_pi_term_exponents(dims_of_all_vars, total_vars, fundamental_dimensions) + pi_terms_data = [] + + # We are considering the repeating variables as starting variables for now. + repeating_idx = [] + for i in 1:fundamental_dimensions + push!(repeating_idx, i) + end + non_repeating_idx = [] # V1 + for i in fundamental_dimensions+1:total_vars + push!(non_repeating_idx, i) + end + + for idx in non_repeating_idx + # Form system of equations for exponents + A = dims_of_all_vars[repeating_idx, 1:fundamental_dimensions]' # Transpose for solving + b = -dims_of_all_vars[idx, 1:fundamental_dimensions] + + exponents = A \ b # Linear solve + + pi_term = Dict("var_idx" => idx, "exponents" => exponents, "repeating_idx" => repeating_idx) + push!(pi_terms_data, pi_term) + end + + return pi_terms_data +end + + +""" + retrieve_pi_terms(Dict{Any}, Vector{Num}) + + Gets the actual PI term provided non repeating index, repeating indices and the original variables +""" +function retrieve_pi_terms(pi_term, var_names) + repeating_idx = pi_term["repeating_idx"] + exp_arr = pi_term["exponents"] + final_pi_term = 1 + + for (ind, val) in enumerate(exp_arr) + exp_arr[ind] = round(val, digits=2) + end + + for (ind, val) in enumerate(repeating_idx) + final_pi_term *= var_names[val]^exp_arr[ind] + end + + # multiplying with non repeating variable for the pi term + final_pi_term *= var_names[pi_term["var_idx"]] + + return final_pi_term +end + + +""" + buckinghumFun(Vector{DynamicQuantities.Quantity}, Vector{Num}) + + Takes an array of DynamicQuantities.Quantity type and variable names separately. Gets the buckinghum + PI terms and returns the array of PI terms +""" +function buckinghumFun(vars_quants, var_names) + + # vars_quants : [m^3/kg, s^-1, ms^-1] + # var_names : [ρ, μ, u, v] + + # Number of variables + total_vars = length(vars_quants) + + # Required for counting fundamental dimensions + # says whicheever units are in picture. In this case : length, mass, time + # [1, 1, 1, 0, 0, 0] + dim_map = zeros(6) + + + # [kg^-3* m, kg*s, s^-1] + dims_vars = [] + for u_arr in vars_quants + push!(dims_vars, u_arr.dimensions) + end + + dims_of_all_vars = get_dims_of_vars(dims_vars, total_vars, dim_map) + + fundamental_dimensions = no_of_fundamental_dims(dim_map) + + pi_terms = find_pi_term_exponents(dims_of_all_vars, total_vars, fundamental_dimensions) + + pis_arr = [] + for val in pi_terms + push!(pis_arr, retrieve_pi_terms(val, var_names)) + end + + return pis_arr +end diff --git a/test/buckinghum.jl b/test/buckinghum.jl new file mode 100644 index 0000000000..dd90ba9fdd --- /dev/null +++ b/test/buckinghum.jl @@ -0,0 +1,21 @@ +include("../src/systems/buckinghum.jl") +using DynamicQuantities +using LinearAlgebra +using Test + +@variables f, d, v, ρ, μ +vars1_arr = [ ρ, d, v, μ, f] + +vars1_quants = [DynamicQuantities.Quantity(0, mass=1, length=-1, time=-1), DynamicQuantities.Quantity(0, length=1) , DynamicQuantities.Quantity(0, length=1, time=-1), DynamicQuantities.Quantity(0, mass=1, length=1, time=-2), DynamicQuantities.Quantity(0, mass=1, length=-1, time=-1)] + +@test isequal(buckinghumFun(vars1_quants, vars1_arr)[1], μ/(d*v*ρ)) +@test isequal(buckinghumFun(vars1_quants, vars1_arr)[2], f/(ρ)) + + +@variables a, b, c, d +vars2_arr = [ a, b , c, d] + +vars2_quants =[DynamicQuantities.Quantity(0, mass=1, length=-3), DynamicQuantities.Quantity(0, mass=1, length=-1, time=-1), DynamicQuantities.Quantity(0, length=1, time=-1), DynamicQuantities.Quantity(0, length=1)] + +@test isequal(buckinghumFun(vars2_quants, vars2_arr)[1], (a*c*d)/b) + From 612f27aae5799a7e9a8281f6e3ce8066ae6577c3 Mon Sep 17 00:00:00 2001 From: Karthik Rishinarada Date: Wed, 5 Mar 2025 12:01:41 +0530 Subject: [PATCH 2/5] style change --- src/systems/buckinghum.jl | 12 +++++------- test/buckinghum.jl | 4 ++-- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/src/systems/buckinghum.jl b/src/systems/buckinghum.jl index e5341e1c1c..52a384eaeb 100644 --- a/src/systems/buckinghum.jl +++ b/src/systems/buckinghum.jl @@ -24,10 +24,10 @@ end get_dims_of_vars(Vector{Any}, Number, Vector{Any}) Gets units of each variable in an array. Returns a matrix where each row corresponds - to binary - Example : kg*m^3 - Each index corresponds to length, mass, time, current, luminosity and temperature - respectively. + to units represented in binary values + Example : kg*m^3 - [3, 1, 0, 0, 0, 0] + In the returned value each row corresponds to length, mass, time, current, luminosity + and temperature respectively. """ function get_dims_of_vars(dims_vars, total_vars, dim_map) # For every single variable it contains row of 1s and 0s mentioning which unit is present @@ -118,9 +118,7 @@ function retrieve_pi_terms(pi_term, var_names) exp_arr = pi_term["exponents"] final_pi_term = 1 - for (ind, val) in enumerate(exp_arr) - exp_arr[ind] = round(val, digits=2) - end + @. exp_arr = round(exp_arr, digits=2) for (ind, val) in enumerate(repeating_idx) final_pi_term *= var_names[val]^exp_arr[ind] diff --git a/test/buckinghum.jl b/test/buckinghum.jl index dd90ba9fdd..716b211a9d 100644 --- a/test/buckinghum.jl +++ b/test/buckinghum.jl @@ -4,7 +4,7 @@ using LinearAlgebra using Test @variables f, d, v, ρ, μ -vars1_arr = [ ρ, d, v, μ, f] +vars1_arr = [ρ, d, v, μ, f] vars1_quants = [DynamicQuantities.Quantity(0, mass=1, length=-1, time=-1), DynamicQuantities.Quantity(0, length=1) , DynamicQuantities.Quantity(0, length=1, time=-1), DynamicQuantities.Quantity(0, mass=1, length=1, time=-2), DynamicQuantities.Quantity(0, mass=1, length=-1, time=-1)] @@ -13,7 +13,7 @@ vars1_quants = [DynamicQuantities.Quantity(0, mass=1, length=-1, time=-1), Dynam @variables a, b, c, d -vars2_arr = [ a, b , c, d] +vars2_arr = [a, b, c, d] vars2_quants =[DynamicQuantities.Quantity(0, mass=1, length=-3), DynamicQuantities.Quantity(0, mass=1, length=-1, time=-1), DynamicQuantities.Quantity(0, length=1, time=-1), DynamicQuantities.Quantity(0, length=1)] From 68f073b69d503e6ba892b2bfffe9b01bf152ad10 Mon Sep 17 00:00:00 2001 From: Karthik Rishinarada Date: Tue, 25 Mar 2025 10:35:08 +0530 Subject: [PATCH 3/5] retrieve new eqs --- src/systems/buckinghum.jl | 26 ++++++++++++++++++++++++++ test/buckinghum.jl | 10 ++++++++-- 2 files changed, 34 insertions(+), 2 deletions(-) diff --git a/src/systems/buckinghum.jl b/src/systems/buckinghum.jl index 52a384eaeb..e96b211364 100644 --- a/src/systems/buckinghum.jl +++ b/src/systems/buckinghum.jl @@ -170,3 +170,29 @@ function buckinghumFun(vars_quants, var_names) return pis_arr end + + +""" + transform_eqs_from_pi_terms(Vector{Any}, Dict{Any, Any}) + + Substitutes the PI terms in the equations to form new set of equations and returns them + in a form of an array. +""" +function transform_eqs_from_pi_terms(pi_eqs, original_equations, iv) + D = ModelingToolkit.Differential(iv) + transformed_eqs = [] + for eq in pi_eqs + pi_term = eq + for each_var in get_variables(pi_term) + println(each_var) + if haskey(original_equations, D(each_var)) != true + original_equations[D(each_var)] = 0 + end + end + + der_eq = expand_derivatives(D(pi_term)) + sub_eq = substitute(der_eq, original_equations) + push!(transformed_eqs, sub_eq) + end + return transformed_eqs +end diff --git a/test/buckinghum.jl b/test/buckinghum.jl index 716b211a9d..804f618de1 100644 --- a/test/buckinghum.jl +++ b/test/buckinghum.jl @@ -12,10 +12,16 @@ vars1_quants = [DynamicQuantities.Quantity(0, mass=1, length=-1, time=-1), Dynam @test isequal(buckinghumFun(vars1_quants, vars1_arr)[2], f/(ρ)) -@variables a, b, c, d +@parameters t, α +@variables a(t), b(t), c(t), d(t) +D = ModelingToolkit.Differential(t) vars2_arr = [a, b, c, d] vars2_quants =[DynamicQuantities.Quantity(0, mass=1, length=-3), DynamicQuantities.Quantity(0, mass=1, length=-1, time=-1), DynamicQuantities.Quantity(0, length=1, time=-1), DynamicQuantities.Quantity(0, length=1)] -@test isequal(buckinghumFun(vars2_quants, vars2_arr)[1], (a*c*d)/b) +pi_terms = buckinghumFun(vars2_quants, vars2_arr) +original_equations = Dict(D(a) => α*t + b, D(b) => c*d/a) + +@test isequal(pi_terms[1], (a*c*d)/b) +@test isequal(transform_eqs_from_pi_terms(pi_terms, original_equations, t)[1], ((b + t*α)*c*d) / b + (-(c^2)*(d^2)) / (b^2)) From aec14ebf12b7bb2d25e0e8073dea711267525489 Mon Sep 17 00:00:00 2001 From: Karthik Rishinarada Date: Tue, 25 Mar 2025 23:48:10 +0530 Subject: [PATCH 4/5] returns new ode system --- src/systems/buckinghum.jl | 32 ++++++++++++++++++++++---------- test/buckinghum.jl | 10 +++++++--- 2 files changed, 29 insertions(+), 13 deletions(-) diff --git a/src/systems/buckinghum.jl b/src/systems/buckinghum.jl index e96b211364..947ddb6442 100644 --- a/src/systems/buckinghum.jl +++ b/src/systems/buckinghum.jl @@ -173,26 +173,38 @@ end """ - transform_eqs_from_pi_terms(Vector{Any}, Dict{Any, Any}) + transform_sys(Vector{Any}, system::ODESystem) - Substitutes the PI terms in the equations to form new set of equations and returns them - in a form of an array. + Substitutes the PI terms in the equations to form new set of equations and returns the + new ODESystem """ -function transform_eqs_from_pi_terms(pi_eqs, original_equations, iv) - D = ModelingToolkit.Differential(iv) + +function transform_sys(pi_eqs, sys::ODESystem, pis_vars) + original_equations = equations(sys) + orginal_eqs_mp = Dict{Any, Any}(eq.lhs => eq.rhs for eq in original_equations) transformed_eqs = [] for eq in pi_eqs pi_term = eq for each_var in get_variables(pi_term) println(each_var) - if haskey(original_equations, D(each_var)) != true - original_equations[D(each_var)] = 0 + if haskey(orginal_eqs_mp, D(each_var)) != true + orginal_eqs_mp[D(each_var)] = 0 end end der_eq = expand_derivatives(D(pi_term)) - sub_eq = substitute(der_eq, original_equations) + sub_eq = substitute(der_eq, orginal_eqs_mp) push!(transformed_eqs, sub_eq) end - return transformed_eqs -end + + equations_for_system = Equation[] + dependent_vars = Any[unknowns(sys)...] + for (i, eq) in pairs(transformed_eqs) + push!(dependent_vars, pis_vars[i]) + push!(equations_for_system, D(pis_vars[i]) ~ eq) + end + + @named new_sys = ODESystem(equations_for_system, ModelingToolkit.get_iv(sys),[π1,a,b,c,d], [α]; defaults=defaults(sys)) + + return new_sys +end \ No newline at end of file diff --git a/test/buckinghum.jl b/test/buckinghum.jl index 804f618de1..a4d56563ca 100644 --- a/test/buckinghum.jl +++ b/test/buckinghum.jl @@ -14,14 +14,18 @@ vars1_quants = [DynamicQuantities.Quantity(0, mass=1, length=-1, time=-1), Dynam @parameters t, α @variables a(t), b(t), c(t), d(t) +@variables π1(t) π2(t) π3(t) π4(t) π5(t) D = ModelingToolkit.Differential(t) vars2_arr = [a, b, c, d] - vars2_quants =[DynamicQuantities.Quantity(0, mass=1, length=-3), DynamicQuantities.Quantity(0, mass=1, length=-1, time=-1), DynamicQuantities.Quantity(0, length=1, time=-1), DynamicQuantities.Quantity(0, length=1)] +original_equations_map = Dict(D(a) => α*t + b, D(b) => c*d/a) +original_equations = [k ~ v for (k,v) in original_equations_map] + +@named sys = ODESystem(original_equations, t, vars2_arr, [α];defaults=Dict(α => 0.5)) pi_terms = buckinghumFun(vars2_quants, vars2_arr) -original_equations = Dict(D(a) => α*t + b, D(b) => c*d/a) +new_sys = transform_sys(pi_terms, sys, [π1, π2, π3, π4, π5]) @test isequal(pi_terms[1], (a*c*d)/b) -@test isequal(transform_eqs_from_pi_terms(pi_terms, original_equations, t)[1], ((b + t*α)*c*d) / b + (-(c^2)*(d^2)) / (b^2)) +@test isequal(equations(new_sys)[1].rhs ,((b + t*α)*c*d) / b + (-(c^2)*(d^2)) / (b^2)) From fc7d2bf50a321381cb723bc6ee06626d107c24dd Mon Sep 17 00:00:00 2001 From: Karthik Rishinarada Date: Tue, 25 Mar 2025 23:49:22 +0530 Subject: [PATCH 5/5] returns new ODESystem after finding pi terms --- src/systems/buckinghum.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/systems/buckinghum.jl b/src/systems/buckinghum.jl index 947ddb6442..8231b127a1 100644 --- a/src/systems/buckinghum.jl +++ b/src/systems/buckinghum.jl @@ -204,7 +204,7 @@ function transform_sys(pi_eqs, sys::ODESystem, pis_vars) push!(equations_for_system, D(pis_vars[i]) ~ eq) end - @named new_sys = ODESystem(equations_for_system, ModelingToolkit.get_iv(sys),[π1,a,b,c,d], [α]; defaults=defaults(sys)) + @named new_sys = ODESystem(equations_for_system, ModelingToolkit.get_iv(sys),dependent_vars, [α]; defaults=defaults(sys)) return new_sys end \ No newline at end of file