BranchesTags

Repository files navigation

TCLB tools

Various tools and extended examples for the TCLB project.

Structure

All tools are divided by language: R, Python, etc. In each directory, the subdirectories should be specific packages related to different aspects of interaction with TCLB

Examples and Papers include scripts used to obtain data used in "real world" cases. Be carefoul - they might be computationally intensive. See Readmi for details

Installation

R

Supervisor: Łukasz Łaniewski-Wołłk

devtools::install_github("CFD-GO/TCLB_tools/R/TCLBtools")

Python

Supervisor: Michał Dzikowski

Content of Python/ directory is meant as a packages repository. Add XXXX/Python to your PYTHONPATH than import packages by subdirectory name.

Python: symbolic_tools

Supervisor: Grzegorz Gruszczyński

It is recomended to use PyCharm to run the scripts. It facilitates recognition of paths and folders. Create a PyCharm project by opening the symbolic_tools/ directory in the editor. Run the examples/ . Some are listed below.

Example: Forcing terms

fromSymbolicCollisions.core.printersimportprint_as_vectorfromSymbolicCollisions.core.ContinousCMTransformsimport \
ContinousCMTransforms, get_mom_vector_from_continuous_deffromSymbolicCollisions.core.cm_symbolsimport \
F3D, dzeta3D, u3DfromSymbolicCollisions.core.cm_symbolsimportrho, moments_dictccmt=ContinousCMTransforms(dzeta3D, u3D, F3D, rho)
lattice='D2Q9'print("\n--- FORCES ---")
print('\n\n// === continuous central moments === \n ')
print('\n//Force -> Force_cm - from continous definition: \n''k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) \n''where fun = forceM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o ')
F_cm=get_mom_vector_from_continuous_def(ccmt.get_force_He_MB,
continuous_transformation=ccmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(F_cm, 'F_cm')

Output:

//Force -> Force_cm - from continous definition: k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) where fun = forceM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o F_cm[0] = 0;
F_cm[1] = Fhydro.x*m00/rho;
F_cm[2] = Fhydro.y*m00/rho;
F_cm[3] = 0;
F_cm[4] = 0;
F_cm[5] = 0;
F_cm[6] = 1/3.*Fhydro.y*m00/rho;
F_cm[7] = 1/3.*Fhydro.x*m00/rho;
F_cm[8] = 0;

Example: Equilibrium Distribution

fromSymbolicCollisions.core.printersimportprint_as_vectorfromsympy.matricesimportMatrixfromsympyimportSymbolfromSymbolicCollisions.core.ContinousCMTransformsimportContinousCMTransforms, get_mom_vector_from_continuous_deffromSymbolicCollisions.core.cm_symbolsimport \
F3D, dzeta3D, u3DfromSymbolicCollisions.core.cm_symbolsimporte_D2Q9, u2D, F2D, rho, moments_dictfromSymbolicCollisions.core.DiscreteCMTransformsimport \
DiscreteCMTransforms, get_mom_vector_from_discrete_deflattice='D2Q9'ccmt=ContinousCMTransforms(dzeta3D, u3D, F3D, rho)
dcmt=DiscreteCMTransforms(e_D2Q9, u2D, F2D, rho)
print("\n--- EQUILIBRIA ---")
print('// === discrete cm ===\n ')
print('\n//population_eq -> cm_eq - by definition: k_mn = sum( (e_ix-ux)^m (e_iy-uy)^n * population_eq_i)')
print("moments: first order (linear) velocity expansion.")
pop_eq=get_mom_vector_from_discrete_def(lambdai: Symbol('m00') *dcmt.get_gamma_first_order(i),
discrete_transform=dcmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(pop_eq, 'cm_eq_first_order')
print("moments: second order (quadratic) velocity expansion.")
pop_eq=get_mom_vector_from_discrete_def(lambdai: Symbol('m00') *dcmt.get_gamma(i),
discrete_transform=dcmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(pop_eq, 'cm_eq_second_order')
print('\n\n// === continous cm === \n ')
# to calculate particular momentrow=moments_dict['D2Q9'][0]
moment=ccmt.get_cm(row, ccmt.get_Maxwellian_DF)
print_as_vector(Matrix([moment]), 'particular_moment')
print('\n//population_eq -> cm_eq - from continous definition: \n''k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) \n''where fun = fM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o ')
cm_eq=get_mom_vector_from_continuous_def(ccmt.get_Maxwellian_DF,
continuous_transformation=ccmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(cm_eq, 'cm_eq')

Output

--- EQUILIBRIA ---
// === discrete cm ===
//population_eq -> cm_eq - by definition: k_mn = sum( (e_ix-ux)^m (e_iy-uy)^n * population_eq_i)
moments: first order (linear) velocity expansion.
cm_eq_first_order[0] = m00;
cm_eq_first_order[1] = 0;
cm_eq_first_order[2] = 0;
cm_eq_first_order[3] = m00*(-ux2 + 1/3.);
cm_eq_first_order[4] = m00*(-uy2 + 1/3.);
cm_eq_first_order[5] = -m00*uxuy;
cm_eq_first_order[6] = 2.*m00*ux2*u.y;
cm_eq_first_order[7] = 2.*m00*u.x*uy2;
cm_eq_first_order[8] = m00*(-3.*ux2*uy2 - 1/3.*ux2 - 1/3.*uy2 + 1/9.);
moments: second order (quadratic) velocity expansion.
cm_eq_second_order[0] = m00;
cm_eq_second_order[1] = 0;
cm_eq_second_order[2] = 0;
cm_eq_second_order[3] = 1/3.*m00;
cm_eq_second_order[4] = 1/3.*m00;
cm_eq_second_order[5] = 0;
cm_eq_second_order[6] = -m00*ux2*u.y;
cm_eq_second_order[7] = -m00*u.x*uy2;
cm_eq_second_order[8] = m00*(3.*ux2*uy2 + 1/9.);
// === continous cm === particular_moment[0] = m00;
//population_eq -> cm_eq - from continous definition: k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) where fun = fM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o cm_eq[0] = m00;
cm_eq[1] = 0;
cm_eq[2] = 0;
cm_eq[3] = 1/3.*m00;
cm_eq[4] = 1/3.*m00;
cm_eq[5] = 0;
cm_eq[6] = 0;
cm_eq[7] = 0;
cm_eq[8] = 1/9.*m00;

Example: Collision Kernel

fromsympy.matricesimporteyefromsympy.printingimportprint_ccodefromSymbolicCollisions.core.cm_symbolsimportomega_ade, omega_b, omega_v, m00fromSymbolicCollisions.core.cm_symbolsimportdynamic_importfromSymbolicCollisions.core.DiscreteCMTransformsimportget_DF, get_m00fromSymbolicCollisions.core.printersimportprint_u2, print_as_vectorfromSymbolicCollisions.core.MatrixGeneratorimportget_raw_moments_matrix, get_shift_matrix# inspired by:# "Consistent Forcing Scheme in the cascaded LBM" L. Fei et al. 2017# eqs 8-12 : (eye(q)-S)*cm + S*cm_eq + (eye(q)-S/2.)*force_in_cm_space# SETUPd=2q=9model='hydro'# choose from '['hydro', 'ade', 'ade_with_f']# DYNAMIC IMPORTSex=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ex_D{d}Q{q}")
ey=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ey_D{d}Q{q}")
ifd==3:
ez=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ez_D{d}Q{q}")
else:
ez=Nonedefget_s_relax_switcher(choice):
s_relax_switcher= {
'hydro': ("SymbolicCollisions.core.cm_symbols", f"S_relax_hydro_D{d}Q{q}"),
'ade_with_f': ("SymbolicCollisions.core.cm_symbols", f"S_relax_ADE_D{d}Q{q}"),
'ade': ("SymbolicCollisions.core.cm_symbols", f"S_relax_ADE_D{d}Q{q}"),
}
which_model=s_relax_switcher.get(choice, lambda: "Invalid argument")
returndynamic_import(*which_model)
S_Relax=get_s_relax_switcher(model)
hardcoded_cm_eq=dynamic_import("SymbolicCollisions.core.hardcoded_results", f"hardcoded_cm_eq_compressible_D{d}Q{q}")
hardcoded_F_cm=dynamic_import("SymbolicCollisions.core.hardcoded_results", f"hardcoded_F_cm_pf_D{d}Q{q}")
fromSymbolicCollisions.core.cm_symbolsimportForce_strasF_str# ARRANGE STUFFMraw=get_raw_moments_matrix(ex, ey, ez)
Nraw=get_shift_matrix(Mraw.inv(), ex, ey, ez)
# from sympy import pprint# pprint(Mraw) # see what you have done# pprint(Nraw)pop_in_str='x_in'# symbol defining populationstemp_pop_str='temp'# symbol defining populationscm_eq_pop_str='cm_eq'# symbol defining populations# GENERATE CODEdefmake_header(choice):
model_switcher= {
'hydro': f"CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t {pop_in_str}[{q}], real_t {omega_v}, vector_t u, vector_t {F_str}) \n{{",
'ade_with_f': f"CudaDeviceFunction void relax_and_collide_ADE_with_F(real_t {pop_in_str}[{q}], real_t {omega_ade}, vector_t u, vector_t {F_str}) \n{{",
'ade': f"CudaDeviceFunction void relax_and_collide_ADE(real_t {pop_in_str}[{q}], real_t {omega_ade}, vector_t u) \n{{",
}
result=model_switcher.get(choice, lambda: "Invalid argument")
print(result)
make_header(model)
print("\t//=== THIS IS AUTOMATICALLY GENERATED CODE ===")
# print(f"real_t {sv} = omega;")# print("real_t bulk_visc = 1./6. ;")# print("real_t {sb} = 1./(3*bulk_visc + 0.5);")# print(f"real_t {sb} = omega_bulk;\n") # s_b = 1./(3*bulk_visc + 0.5)print_u2(d)
print_ccode(get_m00(q, pop_in_str), assign_to=f'\treal_t {m00}')
defmake_variables(choice):
model_switcher= {
'hydro': f"\n\treal_t {temp_pop_str}[{q}];\n",
'ade_with_f': f"\n\treal_t {temp_pop_str}[{q}];\n",
'ade': f"\n\treal_t {temp_pop_str}[{q}];\n",
}
# Get the function from switcher dictionaryresult=model_switcher.get(choice, lambda: "Invalid argument")
print(result)
make_variables(model)
print(f"\tfor (int i = 0; i < {q}; i++) {{\n\t"f"\t{temp_pop_str}[i] = {pop_in_str}[i];}}")
populations=get_DF(q, pop_in_str)
temp_populations=get_DF(q, temp_pop_str)
cm_eq=get_DF(q, cm_eq_pop_str)
F_cm=get_DF(q, F_str)
m=Mraw*temp_populationsprint("\n\t//raw moments from density-probability functions")
# print("\t//[m00, m10, m01, m20, m02, m11, m21, m12, m22]")print_as_vector(m, print_symbol=pop_in_str)
print("\n\t//central moments from raw moments")
cm=Nraw*populationsprint_as_vector(cm, print_symbol=temp_pop_str)
print("\n\t//collision in central moments space")
# print("//calculate equilibrium distributions in cm space")# print("real_t {cm_eq_pop_str}[{q}];\n")# print_as_vector(hardcoded_cm_eq, cm_eq_pop_str) # save time, verbosity# print("//calculate forces in cm space")# print("real_t {F_cm_str}[{q}];")# print_as_vector(hardcoded_F_cm, F_cm_str) # save time, verbosityprint("\t//collide")
defmake_collision(choice):
model_switcher= {
# Relax 2nd moments for hydro, SOI'hydro': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq+ (eye(q) -S_Relax/2) *hardcoded_F_cm,
# Relax 1st moments for ADE, SOI'ade_with_f': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq+ (eye(q) -S_Relax/2) *hardcoded_F_cm,
# Relax 1st moments for ADE, SOI without force'ade': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq,
}
# Get the function from switcher dictionarycm_after_collision=model_switcher.get(choice, lambda: "Invalid argument")
print_as_vector(cm_after_collision, print_symbol=pop_in_str)
make_collision(model)
print("\n\t//back to raw moments")
m=Nraw.inv() *populationsprint_as_vector(m, print_symbol=temp_pop_str)
print("\n\t//back to density-probability functions")
populations=Mraw.inv() *temp_populationsprint_as_vector(populations, print_symbol=pop_in_str)
print("\n}\n")

Output

CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t x_in[9], real_t omega_nu, vector_t u, vector_t F) {
//=== THIS IS AUTOMATICALLY GENERATED CODE ===
real_t uxuy = u.x*u.y;
real_t ux2 = u.x*u.x;
real_t uy2 = u.y*u.y;
real_t m00 = x_in[0] + x_in[1] + x_in[2] + x_in[3] + x_in[4] + x_in[5] + x_in[6] + x_in[7] + x_in[8];
real_t temp[9];for (int i = 0; i < 9; i++) {
temp[i] = x_in[i];}
//raw moments from density-probability functions
x_in[0] = temp[0] + temp[1] + temp[2] + temp[3] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8];
x_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8];
x_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[5] = temp[5] - temp[6] + temp[7] - temp[8];
x_in[6] = temp[5] + temp[6] - temp[7] - temp[8];
x_in[7] = temp[5] - temp[6] - temp[7] + temp[8];
x_in[8] = temp[5] + temp[6] + temp[7] + temp[8];
//central moments from raw moments
temp[0] = x_in[0];
temp[1] = -u.x*x_in[0] + x_in[1];
temp[2] = -u.y*x_in[0] + x_in[2];
temp[3] = ux2*x_in[0] - 2.*u.x*x_in[1] + x_in[3];
temp[4] = uy2*x_in[0] - 2.*u.y*x_in[2] + x_in[4];
temp[5] = uxuy*x_in[0] - u.x*x_in[2] - u.y*x_in[1] + x_in[5];
temp[6] = -ux2*u.y*x_in[0] + ux2*x_in[2] + 2.*uxuy*x_in[1] - 2.*u.x*x_in[5] - u.y*x_in[3] + x_in[6];
temp[7] = -u.x*uy2*x_in[0] + 2.*uxuy*x_in[2] - u.x*x_in[4] + uy2*x_in[1] - 2.*u.y*x_in[5] + x_in[7];
temp[8] = ux2*uy2*x_in[0] - 2.*ux2*u.y*x_in[2] + ux2*x_in[4] - 2.*u.x*uy2*x_in[1] + 4.*uxuy*x_in[5] - 2.*u.x*x_in[7] + uy2*x_in[3] - 2.*u.y*x_in[6] + x_in[8];
//collision in central moments space
//collide
x_in[0] = m00;
x_in[1] = 1/2.*F.x;
x_in[2] = 1/2.*F.y;
x_in[3] = 1/3.*m00*omega_bulk - 1/2.*omega_bulk*temp[3] - 1/2.*omega_bulk*temp[4] - 1/2.*omega_nu*temp[3] + 1/2.*omega_nu*temp[4] + temp[3];
x_in[4] = 1/3.*m00*omega_bulk - 1/2.*omega_bulk*temp[3] - 1/2.*omega_bulk*temp[4] + 1/2.*omega_nu*temp[3] - 1/2.*omega_nu*temp[4] + temp[4];
x_in[5] = -temp[5]*(omega_nu - 1.);
x_in[6] = 1/6.*F.y;
x_in[7] = 1/6.*F.x;
x_in[8] = 1/9.*m00;
//back to raw moments
temp[0] = x_in[0];
temp[1] = u.x*x_in[0] + x_in[1];
temp[2] = u.y*x_in[0] + x_in[2];
temp[3] = ux2*x_in[0] + 2.*u.x*x_in[1] + x_in[3];
temp[4] = uy2*x_in[0] + 2.*u.y*x_in[2] + x_in[4];
temp[5] = uxuy*x_in[0] + u.x*x_in[2] + u.y*x_in[1] + x_in[5];
temp[6] = ux2*u.y*x_in[0] + ux2*x_in[2] + 2.*uxuy*x_in[1] + 2.*u.x*x_in[5] + u.y*x_in[3] + x_in[6];
temp[7] = u.x*uy2*x_in[0] + 2.*uxuy*x_in[2] + u.x*x_in[4] + uy2*x_in[1] + 2.*u.y*x_in[5] + x_in[7];
temp[8] = ux2*uy2*x_in[0] + 2.*ux2*u.y*x_in[2] + ux2*x_in[4] + 2.*u.x*uy2*x_in[1] + 4.*uxuy*x_in[5] + 2.*u.x*x_in[7] + uy2*x_in[3] + 2.*u.y*x_in[6] + x_in[8];
//back to density-probability functions
x_in[0] = temp[0] - temp[3] - temp[4] + temp[8];
x_in[1] = 1/2.*temp[1] + 1/2.*temp[3] - 1/2.*temp[7] - 1/2.*temp[8];
x_in[2] = 1/2.*temp[2] + 1/2.*temp[4] - 1/2.*temp[6] - 1/2.*temp[8];
x_in[3] = -1/2.*temp[1] + 1/2.*temp[3] + 1/2.*temp[7] - 1/2.*temp[8];
x_in[4] = -1/2.*temp[2] + 1/2.*temp[4] + 1/2.*temp[6] - 1/2.*temp[8];
x_in[5] = 1/4.*temp[5] + 1/4.*temp[6] + 1/4.*temp[7] + 1/4.*temp[8];
x_in[6] = -1/4.*temp[5] + 1/4.*temp[6] - 1/4.*temp[7] + 1/4.*temp[8];
x_in[7] = 1/4.*temp[5] - 1/4.*temp[6] - 1/4.*temp[7] + 1/4.*temp[8];
x_in[8] = -1/4.*temp[5] - 1/4.*temp[6] + 1/4.*temp[7] + 1/4.*temp[8];
}

About

Various tools for the TCLB project

Resources

Readme
GPL-3.0 license
Activity

Stars

0 stars

Watchers

1 watching

Forks

0 forks
Report repository

Releases

Packages

Contributors

Languages

, 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Add copy buttons to all 
 blocks\n(function() {\n function addCopyButtons() {\n document.querySelectorAll('pre code').forEach(function(codeBlock) {\n if (codeBlock.parentElement.hasAttribute('data-copy-added')) return;\n codeBlock.parentElement.setAttribute('data-copy-added', 'true');\n \n var btn = document.createElement('button');\n btn.textContent = 'Copy';\n btn.style.cssText = 'position:absolute;top:4px;right:4px;padding:2px 8px;font-size:11px;background:#4ecdc4;border:none;border-radius:4px;color:#1a1a2e;cursor:pointer;opacity:0.7;transition:opacity 0.2s;';\n btn.onmouseover = function() { this.style.opacity = '1'; };\n btn.onmouseout = function() { this.style.opacity = '0.7'; };\n btn.onclick = function() {\n navigator.clipboard.writeText(codeBlock.textContent).then(function() {\n btn.textContent = 'Copied!';\n setTimeout(function() { btn.textContent = 'Copy'; }, 1500);\n });\n };\n codeBlock.parentElement.style.position = 'relative';\n codeBlock.parentElement.appendChild(btn);\n });\n }\n \n addCopyButtons();\n \n // Re-run on dynamic content\n var observer = new MutationObserver(addCopyButtons);\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "Add Copy Buttons to Code Blocks");
}
} catch(__e) { console.warn('[Userscript:Add Copy Buttons to Code Blocks]', __e); }
})();
(function(){
try {
var __m = "github.com";
var __re = new RegExp('^' + "github\\.com" + '

Skip to content
BranchesTags
Repository files navigation
TCLB tools
Various tools and extended examples for the TCLB project.
Structure
All tools are divided by language: R, Python, etc. In each directory, the subdirectories should be specific packages related to different aspects of interaction with TCLB
Examples and Papers include scripts used to obtain data used in "real world" cases. Be carefoul - they might be computationally intensive. See Readmi for details
Installation
R
Supervisor: Łukasz Łaniewski-Wołłk
devtools::install_github("CFD-GO/TCLB_tools/R/TCLBtools")
Python
Supervisor: Michał Dzikowski
Content of Python/ directory is meant as a packages repository. Add XXXX/Python to your PYTHONPATH than import packages by subdirectory name.
Python: symbolic_tools
Supervisor: Grzegorz Gruszczyński
It is recomended to use PyCharm to run the scripts.
It facilitates recognition of paths and folders.
Create a PyCharm project by opening the symbolic_tools/ directory in the editor.
Run the examples/ . Some are listed below.
Example: Forcing terms
fromSymbolicCollisions.core.printersimportprint_as_vectorfromSymbolicCollisions.core.ContinousCMTransformsimport \
ContinousCMTransforms, get_mom_vector_from_continuous_deffromSymbolicCollisions.core.cm_symbolsimport \
F3D, dzeta3D, u3DfromSymbolicCollisions.core.cm_symbolsimportrho, moments_dictccmt=ContinousCMTransforms(dzeta3D, u3D, F3D, rho)
lattice='D2Q9'print("\n--- FORCES ---")
print('\n\n// === continuous central moments === \n ')
print('\n//Force -> Force_cm - from continous definition: \n''k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) \n''where fun = forceM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o ')
F_cm=get_mom_vector_from_continuous_def(ccmt.get_force_He_MB,
continuous_transformation=ccmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(F_cm, 'F_cm')
Output:
//Force -> Force_cm - from continous definition: k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) where fun = forceM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o F_cm[0] = 0;
F_cm[1] = Fhydro.x*m00/rho;
F_cm[2] = Fhydro.y*m00/rho;
F_cm[3] = 0;
F_cm[4] = 0;
F_cm[5] = 0;
F_cm[6] = 1/3.*Fhydro.y*m00/rho;
F_cm[7] = 1/3.*Fhydro.x*m00/rho;
F_cm[8] = 0;
Example: Equilibrium Distribution
fromSymbolicCollisions.core.printersimportprint_as_vectorfromsympy.matricesimportMatrixfromsympyimportSymbolfromSymbolicCollisions.core.ContinousCMTransformsimportContinousCMTransforms, get_mom_vector_from_continuous_deffromSymbolicCollisions.core.cm_symbolsimport \
F3D, dzeta3D, u3DfromSymbolicCollisions.core.cm_symbolsimporte_D2Q9, u2D, F2D, rho, moments_dictfromSymbolicCollisions.core.DiscreteCMTransformsimport \
DiscreteCMTransforms, get_mom_vector_from_discrete_deflattice='D2Q9'ccmt=ContinousCMTransforms(dzeta3D, u3D, F3D, rho)
dcmt=DiscreteCMTransforms(e_D2Q9, u2D, F2D, rho)
print("\n--- EQUILIBRIA ---")
print('// === discrete cm ===\n ')
print('\n//population_eq -> cm_eq - by definition: k_mn = sum( (e_ix-ux)^m (e_iy-uy)^n * population_eq_i)')
print("moments: first order (linear) velocity expansion.")
pop_eq=get_mom_vector_from_discrete_def(lambdai: Symbol('m00') *dcmt.get_gamma_first_order(i),
discrete_transform=dcmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(pop_eq, 'cm_eq_first_order')
print("moments: second order (quadratic) velocity expansion.")
pop_eq=get_mom_vector_from_discrete_def(lambdai: Symbol('m00') *dcmt.get_gamma(i),
discrete_transform=dcmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(pop_eq, 'cm_eq_second_order')
print('\n\n// === continous cm === \n ')
# to calculate particular momentrow=moments_dict['D2Q9'][0]
moment=ccmt.get_cm(row, ccmt.get_Maxwellian_DF)
print_as_vector(Matrix([moment]), 'particular_moment')
print('\n//population_eq -> cm_eq - from continous definition: \n''k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) \n''where fun = fM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o ')
cm_eq=get_mom_vector_from_continuous_def(ccmt.get_Maxwellian_DF,
continuous_transformation=ccmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(cm_eq, 'cm_eq')
Output
--- EQUILIBRIA ---
// === discrete cm ===
//population_eq -> cm_eq - by definition: k_mn = sum( (e_ix-ux)^m (e_iy-uy)^n * population_eq_i)
moments: first order (linear) velocity expansion.
cm_eq_first_order[0] = m00;
cm_eq_first_order[1] = 0;
cm_eq_first_order[2] = 0;
cm_eq_first_order[3] = m00*(-ux2 + 1/3.);
cm_eq_first_order[4] = m00*(-uy2 + 1/3.);
cm_eq_first_order[5] = -m00*uxuy;
cm_eq_first_order[6] = 2.*m00*ux2*u.y;
cm_eq_first_order[7] = 2.*m00*u.x*uy2;
cm_eq_first_order[8] = m00*(-3.*ux2*uy2 - 1/3.*ux2 - 1/3.*uy2 + 1/9.);
moments: second order (quadratic) velocity expansion.
cm_eq_second_order[0] = m00;
cm_eq_second_order[1] = 0;
cm_eq_second_order[2] = 0;
cm_eq_second_order[3] = 1/3.*m00;
cm_eq_second_order[4] = 1/3.*m00;
cm_eq_second_order[5] = 0;
cm_eq_second_order[6] = -m00*ux2*u.y;
cm_eq_second_order[7] = -m00*u.x*uy2;
cm_eq_second_order[8] = m00*(3.*ux2*uy2 + 1/9.);
// === continous cm === particular_moment[0] = m00;
//population_eq -> cm_eq - from continous definition: k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) where fun = fM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o cm_eq[0] = m00;
cm_eq[1] = 0;
cm_eq[2] = 0;
cm_eq[3] = 1/3.*m00;
cm_eq[4] = 1/3.*m00;
cm_eq[5] = 0;
cm_eq[6] = 0;
cm_eq[7] = 0;
cm_eq[8] = 1/9.*m00;
Example: Collision Kernel
fromsympy.matricesimporteyefromsympy.printingimportprint_ccodefromSymbolicCollisions.core.cm_symbolsimportomega_ade, omega_b, omega_v, m00fromSymbolicCollisions.core.cm_symbolsimportdynamic_importfromSymbolicCollisions.core.DiscreteCMTransformsimportget_DF, get_m00fromSymbolicCollisions.core.printersimportprint_u2, print_as_vectorfromSymbolicCollisions.core.MatrixGeneratorimportget_raw_moments_matrix, get_shift_matrix# inspired by:# "Consistent Forcing Scheme in the cascaded LBM" L. Fei et al. 2017# eqs 8-12 : (eye(q)-S)*cm + S*cm_eq + (eye(q)-S/2.)*force_in_cm_space# SETUPd=2q=9model='hydro'# choose from '['hydro', 'ade', 'ade_with_f']# DYNAMIC IMPORTSex=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ex_D{d}Q{q}")
ey=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ey_D{d}Q{q}")
ifd==3:
ez=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ez_D{d}Q{q}")
else:
ez=Nonedefget_s_relax_switcher(choice):
s_relax_switcher= {
'hydro': ("SymbolicCollisions.core.cm_symbols", f"S_relax_hydro_D{d}Q{q}"),
'ade_with_f': ("SymbolicCollisions.core.cm_symbols", f"S_relax_ADE_D{d}Q{q}"),
'ade': ("SymbolicCollisions.core.cm_symbols", f"S_relax_ADE_D{d}Q{q}"),
}
which_model=s_relax_switcher.get(choice, lambda: "Invalid argument")
returndynamic_import(*which_model)
S_Relax=get_s_relax_switcher(model)
hardcoded_cm_eq=dynamic_import("SymbolicCollisions.core.hardcoded_results", f"hardcoded_cm_eq_compressible_D{d}Q{q}")
hardcoded_F_cm=dynamic_import("SymbolicCollisions.core.hardcoded_results", f"hardcoded_F_cm_pf_D{d}Q{q}")
fromSymbolicCollisions.core.cm_symbolsimportForce_strasF_str# ARRANGE STUFFMraw=get_raw_moments_matrix(ex, ey, ez)
Nraw=get_shift_matrix(Mraw.inv(), ex, ey, ez)
# from sympy import pprint# pprint(Mraw) # see what you have done# pprint(Nraw)pop_in_str='x_in'# symbol defining populationstemp_pop_str='temp'# symbol defining populationscm_eq_pop_str='cm_eq'# symbol defining populations# GENERATE CODEdefmake_header(choice):
model_switcher= {
'hydro': f"CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t {pop_in_str}[{q}], real_t {omega_v}, vector_t u, vector_t {F_str}) \n{{",
'ade_with_f': f"CudaDeviceFunction void relax_and_collide_ADE_with_F(real_t {pop_in_str}[{q}], real_t {omega_ade}, vector_t u, vector_t {F_str}) \n{{",
'ade': f"CudaDeviceFunction void relax_and_collide_ADE(real_t {pop_in_str}[{q}], real_t {omega_ade}, vector_t u) \n{{",
}
result=model_switcher.get(choice, lambda: "Invalid argument")
print(result)
make_header(model)
print("\t//=== THIS IS AUTOMATICALLY GENERATED CODE ===")
# print(f"real_t {sv} = omega;")# print("real_t bulk_visc = 1./6. ;")# print("real_t {sb} = 1./(3*bulk_visc + 0.5);")# print(f"real_t {sb} = omega_bulk;\n") # s_b = 1./(3*bulk_visc + 0.5)print_u2(d)
print_ccode(get_m00(q, pop_in_str), assign_to=f'\treal_t {m00}')
defmake_variables(choice):
model_switcher= {
'hydro': f"\n\treal_t {temp_pop_str}[{q}];\n",
'ade_with_f': f"\n\treal_t {temp_pop_str}[{q}];\n",
'ade': f"\n\treal_t {temp_pop_str}[{q}];\n",
}
# Get the function from switcher dictionaryresult=model_switcher.get(choice, lambda: "Invalid argument")
print(result)
make_variables(model)
print(f"\tfor (int i = 0; i < {q}; i++) {{\n\t"f"\t{temp_pop_str}[i] = {pop_in_str}[i];}}")
populations=get_DF(q, pop_in_str)
temp_populations=get_DF(q, temp_pop_str)
cm_eq=get_DF(q, cm_eq_pop_str)
F_cm=get_DF(q, F_str)
m=Mraw*temp_populationsprint("\n\t//raw moments from density-probability functions")
# print("\t//[m00, m10, m01, m20, m02, m11, m21, m12, m22]")print_as_vector(m, print_symbol=pop_in_str)
print("\n\t//central moments from raw moments")
cm=Nraw*populationsprint_as_vector(cm, print_symbol=temp_pop_str)
print("\n\t//collision in central moments space")
# print("//calculate equilibrium distributions in cm space")# print("real_t {cm_eq_pop_str}[{q}];\n")# print_as_vector(hardcoded_cm_eq, cm_eq_pop_str) # save time, verbosity# print("//calculate forces in cm space")# print("real_t {F_cm_str}[{q}];")# print_as_vector(hardcoded_F_cm, F_cm_str) # save time, verbosityprint("\t//collide")
defmake_collision(choice):
model_switcher= {
# Relax 2nd moments for hydro, SOI'hydro': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq+ (eye(q) -S_Relax/2) *hardcoded_F_cm,
# Relax 1st moments for ADE, SOI'ade_with_f': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq+ (eye(q) -S_Relax/2) *hardcoded_F_cm,
# Relax 1st moments for ADE, SOI without force'ade': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq,
}
# Get the function from switcher dictionarycm_after_collision=model_switcher.get(choice, lambda: "Invalid argument")
print_as_vector(cm_after_collision, print_symbol=pop_in_str)
make_collision(model)
print("\n\t//back to raw moments")
m=Nraw.inv() *populationsprint_as_vector(m, print_symbol=temp_pop_str)
print("\n\t//back to density-probability functions")
populations=Mraw.inv() *temp_populationsprint_as_vector(populations, print_symbol=pop_in_str)
print("\n}\n")
Output
CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t x_in[9], real_t omega_nu, vector_t u, vector_t F) {
//=== THIS IS AUTOMATICALLY GENERATED CODE ===
real_t uxuy = u.x*u.y;
real_t ux2 = u.x*u.x;
real_t uy2 = u.y*u.y;
real_t m00 = x_in[0] + x_in[1] + x_in[2] + x_in[3] + x_in[4] + x_in[5] + x_in[6] + x_in[7] + x_in[8];
real_t temp[9];for (int i = 0; i < 9; i++) {
temp[i] = x_in[i];}
//raw moments from density-probability functions
x_in[0] = temp[0] + temp[1] + temp[2] + temp[3] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8];
x_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8];
x_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[5] = temp[5] - temp[6] + temp[7] - temp[8];
x_in[6] = temp[5] + temp[6] - temp[7] - temp[8];
x_in[7] = temp[5] - temp[6] - temp[7] + temp[8];
x_in[8] = temp[5] + temp[6] + temp[7] + temp[8];
//central moments from raw moments
temp[0] = x_in[0];
temp[1] = -u.x*x_in[0] + x_in[1];
temp[2] = -u.y*x_in[0] + x_in[2];
temp[3] = ux2*x_in[0] - 2.*u.x*x_in[1] + x_in[3];
temp[4] = uy2*x_in[0] - 2.*u.y*x_in[2] + x_in[4];
temp[5] = uxuy*x_in[0] - u.x*x_in[2] - u.y*x_in[1] + x_in[5];
temp[6] = -ux2*u.y*x_in[0] + ux2*x_in[2] + 2.*uxuy*x_in[1] - 2.*u.x*x_in[5] - u.y*x_in[3] + x_in[6];
temp[7] = -u.x*uy2*x_in[0] + 2.*uxuy*x_in[2] - u.x*x_in[4] + uy2*x_in[1] - 2.*u.y*x_in[5] + x_in[7];
temp[8] = ux2*uy2*x_in[0] - 2.*ux2*u.y*x_in[2] + ux2*x_in[4] - 2.*u.x*uy2*x_in[1] + 4.*uxuy*x_in[5] - 2.*u.x*x_in[7] + uy2*x_in[3] - 2.*u.y*x_in[6] + x_in[8];
//collision in central moments space
//collide
x_in[0] = m00;
x_in[1] = 1/2.*F.x;
x_in[2] = 1/2.*F.y;
x_in[3] = 1/3.*m00*omega_bulk - 1/2.*omega_bulk*temp[3] - 1/2.*omega_bulk*temp[4] - 1/2.*omega_nu*temp[3] + 1/2.*omega_nu*temp[4] + temp[3];
x_in[4] = 1/3.*m00*omega_bulk - 1/2.*omega_bulk*temp[3] - 1/2.*omega_bulk*temp[4] + 1/2.*omega_nu*temp[3] - 1/2.*omega_nu*temp[4] + temp[4];
x_in[5] = -temp[5]*(omega_nu - 1.);
x_in[6] = 1/6.*F.y;
x_in[7] = 1/6.*F.x;
x_in[8] = 1/9.*m00;
//back to raw moments
temp[0] = x_in[0];
temp[1] = u.x*x_in[0] + x_in[1];
temp[2] = u.y*x_in[0] + x_in[2];
temp[3] = ux2*x_in[0] + 2.*u.x*x_in[1] + x_in[3];
temp[4] = uy2*x_in[0] + 2.*u.y*x_in[2] + x_in[4];
temp[5] = uxuy*x_in[0] + u.x*x_in[2] + u.y*x_in[1] + x_in[5];
temp[6] = ux2*u.y*x_in[0] + ux2*x_in[2] + 2.*uxuy*x_in[1] + 2.*u.x*x_in[5] + u.y*x_in[3] + x_in[6];
temp[7] = u.x*uy2*x_in[0] + 2.*uxuy*x_in[2] + u.x*x_in[4] + uy2*x_in[1] + 2.*u.y*x_in[5] + x_in[7];
temp[8] = ux2*uy2*x_in[0] + 2.*ux2*u.y*x_in[2] + ux2*x_in[4] + 2.*u.x*uy2*x_in[1] + 4.*uxuy*x_in[5] + 2.*u.x*x_in[7] + uy2*x_in[3] + 2.*u.y*x_in[6] + x_in[8];
//back to density-probability functions
x_in[0] = temp[0] - temp[3] - temp[4] + temp[8];
x_in[1] = 1/2.*temp[1] + 1/2.*temp[3] - 1/2.*temp[7] - 1/2.*temp[8];
x_in[2] = 1/2.*temp[2] + 1/2.*temp[4] - 1/2.*temp[6] - 1/2.*temp[8];
x_in[3] = -1/2.*temp[1] + 1/2.*temp[3] + 1/2.*temp[7] - 1/2.*temp[8];
x_in[4] = -1/2.*temp[2] + 1/2.*temp[4] + 1/2.*temp[6] - 1/2.*temp[8];
x_in[5] = 1/4.*temp[5] + 1/4.*temp[6] + 1/4.*temp[7] + 1/4.*temp[8];
x_in[6] = -1/4.*temp[5] + 1/4.*temp[6] - 1/4.*temp[7] + 1/4.*temp[8];
x_in[7] = 1/4.*temp[5] - 1/4.*temp[6] - 1/4.*temp[7] + 1/4.*temp[8];
x_in[8] = -1/4.*temp[5] - 1/4.*temp[6] + 1/4.*temp[7] + 1/4.*temp[8];
}
About
Various tools for the TCLB project
Resources
Readme
GPL-3.0 license
Activity
Stars
0 stars
Watchers
1 watching
Forks
0 forks
Report repository
Releases
Packages
Contributors
Languages

, 'i');
if (__m === '*' || __re.test(location.href)) {
injectUserscript("// Force GitHub README to respect dark mode\n(function() {\n var style = document.createElement('style');\n style.textContent = '\n .markdown-body {\n color-scheme: dark light;\n }\n .markdown-body pre { background: #161b22 !important; }\n .markdown-body code { background: rgba(110, 118, 129, 0.4) !important; }\n .markdown-body table th, .markdown-body table td { border-color: #30363d !important; }\n .markdown-body img { background: #0d1117; }\n .markdown-body blockquote { border-left-color: #8b949e; }\n .markdown-body hr { border-color: #30363d; }\n ';\n document.head.appendChild(style);\n})();", "GitHub Dark Mode README Fix");
}
} catch(__e) { console.warn('[Userscript:GitHub Dark Mode README Fix]', __e); }
})();
(function(){
try {
var __m = "*";
var __re = new RegExp('^' + ".*" + '

Skip to content
BranchesTags
Repository files navigation
TCLB tools
Various tools and extended examples for the TCLB project.
Structure
All tools are divided by language: R, Python, etc. In each directory, the subdirectories should be specific packages related to different aspects of interaction with TCLB
Examples and Papers include scripts used to obtain data used in "real world" cases. Be carefoul - they might be computationally intensive. See Readmi for details
Installation
R
Supervisor: Łukasz Łaniewski-Wołłk
devtools::install_github("CFD-GO/TCLB_tools/R/TCLBtools")
Python
Supervisor: Michał Dzikowski
Content of Python/ directory is meant as a packages repository. Add XXXX/Python to your PYTHONPATH than import packages by subdirectory name.
Python: symbolic_tools
Supervisor: Grzegorz Gruszczyński
It is recomended to use PyCharm to run the scripts.
It facilitates recognition of paths and folders.
Create a PyCharm project by opening the symbolic_tools/ directory in the editor.
Run the examples/ . Some are listed below.
Example: Forcing terms
fromSymbolicCollisions.core.printersimportprint_as_vectorfromSymbolicCollisions.core.ContinousCMTransformsimport \
ContinousCMTransforms, get_mom_vector_from_continuous_deffromSymbolicCollisions.core.cm_symbolsimport \
F3D, dzeta3D, u3DfromSymbolicCollisions.core.cm_symbolsimportrho, moments_dictccmt=ContinousCMTransforms(dzeta3D, u3D, F3D, rho)
lattice='D2Q9'print("\n--- FORCES ---")
print('\n\n// === continuous central moments === \n ')
print('\n//Force -> Force_cm - from continous definition: \n''k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) \n''where fun = forceM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o ')
F_cm=get_mom_vector_from_continuous_def(ccmt.get_force_He_MB,
continuous_transformation=ccmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(F_cm, 'F_cm')
Output:
//Force -> Force_cm - from continous definition: k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) where fun = forceM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o F_cm[0] = 0;
F_cm[1] = Fhydro.x*m00/rho;
F_cm[2] = Fhydro.y*m00/rho;
F_cm[3] = 0;
F_cm[4] = 0;
F_cm[5] = 0;
F_cm[6] = 1/3.*Fhydro.y*m00/rho;
F_cm[7] = 1/3.*Fhydro.x*m00/rho;
F_cm[8] = 0;
Example: Equilibrium Distribution
fromSymbolicCollisions.core.printersimportprint_as_vectorfromsympy.matricesimportMatrixfromsympyimportSymbolfromSymbolicCollisions.core.ContinousCMTransformsimportContinousCMTransforms, get_mom_vector_from_continuous_deffromSymbolicCollisions.core.cm_symbolsimport \
F3D, dzeta3D, u3DfromSymbolicCollisions.core.cm_symbolsimporte_D2Q9, u2D, F2D, rho, moments_dictfromSymbolicCollisions.core.DiscreteCMTransformsimport \
DiscreteCMTransforms, get_mom_vector_from_discrete_deflattice='D2Q9'ccmt=ContinousCMTransforms(dzeta3D, u3D, F3D, rho)
dcmt=DiscreteCMTransforms(e_D2Q9, u2D, F2D, rho)
print("\n--- EQUILIBRIA ---")
print('// === discrete cm ===\n ')
print('\n//population_eq -> cm_eq - by definition: k_mn = sum( (e_ix-ux)^m (e_iy-uy)^n * population_eq_i)')
print("moments: first order (linear) velocity expansion.")
pop_eq=get_mom_vector_from_discrete_def(lambdai: Symbol('m00') *dcmt.get_gamma_first_order(i),
discrete_transform=dcmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(pop_eq, 'cm_eq_first_order')
print("moments: second order (quadratic) velocity expansion.")
pop_eq=get_mom_vector_from_discrete_def(lambdai: Symbol('m00') *dcmt.get_gamma(i),
discrete_transform=dcmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(pop_eq, 'cm_eq_second_order')
print('\n\n// === continous cm === \n ')
# to calculate particular momentrow=moments_dict['D2Q9'][0]
moment=ccmt.get_cm(row, ccmt.get_Maxwellian_DF)
print_as_vector(Matrix([moment]), 'particular_moment')
print('\n//population_eq -> cm_eq - from continous definition: \n''k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) \n''where fun = fM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o ')
cm_eq=get_mom_vector_from_continuous_def(ccmt.get_Maxwellian_DF,
continuous_transformation=ccmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(cm_eq, 'cm_eq')
Output
--- EQUILIBRIA ---
// === discrete cm ===
//population_eq -> cm_eq - by definition: k_mn = sum( (e_ix-ux)^m (e_iy-uy)^n * population_eq_i)
moments: first order (linear) velocity expansion.
cm_eq_first_order[0] = m00;
cm_eq_first_order[1] = 0;
cm_eq_first_order[2] = 0;
cm_eq_first_order[3] = m00*(-ux2 + 1/3.);
cm_eq_first_order[4] = m00*(-uy2 + 1/3.);
cm_eq_first_order[5] = -m00*uxuy;
cm_eq_first_order[6] = 2.*m00*ux2*u.y;
cm_eq_first_order[7] = 2.*m00*u.x*uy2;
cm_eq_first_order[8] = m00*(-3.*ux2*uy2 - 1/3.*ux2 - 1/3.*uy2 + 1/9.);
moments: second order (quadratic) velocity expansion.
cm_eq_second_order[0] = m00;
cm_eq_second_order[1] = 0;
cm_eq_second_order[2] = 0;
cm_eq_second_order[3] = 1/3.*m00;
cm_eq_second_order[4] = 1/3.*m00;
cm_eq_second_order[5] = 0;
cm_eq_second_order[6] = -m00*ux2*u.y;
cm_eq_second_order[7] = -m00*u.x*uy2;
cm_eq_second_order[8] = m00*(3.*ux2*uy2 + 1/9.);
// === continous cm === particular_moment[0] = m00;
//population_eq -> cm_eq - from continous definition: k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) where fun = fM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o cm_eq[0] = m00;
cm_eq[1] = 0;
cm_eq[2] = 0;
cm_eq[3] = 1/3.*m00;
cm_eq[4] = 1/3.*m00;
cm_eq[5] = 0;
cm_eq[6] = 0;
cm_eq[7] = 0;
cm_eq[8] = 1/9.*m00;
Example: Collision Kernel
fromsympy.matricesimporteyefromsympy.printingimportprint_ccodefromSymbolicCollisions.core.cm_symbolsimportomega_ade, omega_b, omega_v, m00fromSymbolicCollisions.core.cm_symbolsimportdynamic_importfromSymbolicCollisions.core.DiscreteCMTransformsimportget_DF, get_m00fromSymbolicCollisions.core.printersimportprint_u2, print_as_vectorfromSymbolicCollisions.core.MatrixGeneratorimportget_raw_moments_matrix, get_shift_matrix# inspired by:# "Consistent Forcing Scheme in the cascaded LBM" L. Fei et al. 2017# eqs 8-12 : (eye(q)-S)*cm + S*cm_eq + (eye(q)-S/2.)*force_in_cm_space# SETUPd=2q=9model='hydro'# choose from '['hydro', 'ade', 'ade_with_f']# DYNAMIC IMPORTSex=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ex_D{d}Q{q}")
ey=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ey_D{d}Q{q}")
ifd==3:
ez=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ez_D{d}Q{q}")
else:
ez=Nonedefget_s_relax_switcher(choice):
s_relax_switcher= {
'hydro': ("SymbolicCollisions.core.cm_symbols", f"S_relax_hydro_D{d}Q{q}"),
'ade_with_f': ("SymbolicCollisions.core.cm_symbols", f"S_relax_ADE_D{d}Q{q}"),
'ade': ("SymbolicCollisions.core.cm_symbols", f"S_relax_ADE_D{d}Q{q}"),
}
which_model=s_relax_switcher.get(choice, lambda: "Invalid argument")
returndynamic_import(*which_model)
S_Relax=get_s_relax_switcher(model)
hardcoded_cm_eq=dynamic_import("SymbolicCollisions.core.hardcoded_results", f"hardcoded_cm_eq_compressible_D{d}Q{q}")
hardcoded_F_cm=dynamic_import("SymbolicCollisions.core.hardcoded_results", f"hardcoded_F_cm_pf_D{d}Q{q}")
fromSymbolicCollisions.core.cm_symbolsimportForce_strasF_str# ARRANGE STUFFMraw=get_raw_moments_matrix(ex, ey, ez)
Nraw=get_shift_matrix(Mraw.inv(), ex, ey, ez)
# from sympy import pprint# pprint(Mraw) # see what you have done# pprint(Nraw)pop_in_str='x_in'# symbol defining populationstemp_pop_str='temp'# symbol defining populationscm_eq_pop_str='cm_eq'# symbol defining populations# GENERATE CODEdefmake_header(choice):
model_switcher= {
'hydro': f"CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t {pop_in_str}[{q}], real_t {omega_v}, vector_t u, vector_t {F_str}) \n{{",
'ade_with_f': f"CudaDeviceFunction void relax_and_collide_ADE_with_F(real_t {pop_in_str}[{q}], real_t {omega_ade}, vector_t u, vector_t {F_str}) \n{{",
'ade': f"CudaDeviceFunction void relax_and_collide_ADE(real_t {pop_in_str}[{q}], real_t {omega_ade}, vector_t u) \n{{",
}
result=model_switcher.get(choice, lambda: "Invalid argument")
print(result)
make_header(model)
print("\t//=== THIS IS AUTOMATICALLY GENERATED CODE ===")
# print(f"real_t {sv} = omega;")# print("real_t bulk_visc = 1./6. ;")# print("real_t {sb} = 1./(3*bulk_visc + 0.5);")# print(f"real_t {sb} = omega_bulk;\n") # s_b = 1./(3*bulk_visc + 0.5)print_u2(d)
print_ccode(get_m00(q, pop_in_str), assign_to=f'\treal_t {m00}')
defmake_variables(choice):
model_switcher= {
'hydro': f"\n\treal_t {temp_pop_str}[{q}];\n",
'ade_with_f': f"\n\treal_t {temp_pop_str}[{q}];\n",
'ade': f"\n\treal_t {temp_pop_str}[{q}];\n",
}
# Get the function from switcher dictionaryresult=model_switcher.get(choice, lambda: "Invalid argument")
print(result)
make_variables(model)
print(f"\tfor (int i = 0; i < {q}; i++) {{\n\t"f"\t{temp_pop_str}[i] = {pop_in_str}[i];}}")
populations=get_DF(q, pop_in_str)
temp_populations=get_DF(q, temp_pop_str)
cm_eq=get_DF(q, cm_eq_pop_str)
F_cm=get_DF(q, F_str)
m=Mraw*temp_populationsprint("\n\t//raw moments from density-probability functions")
# print("\t//[m00, m10, m01, m20, m02, m11, m21, m12, m22]")print_as_vector(m, print_symbol=pop_in_str)
print("\n\t//central moments from raw moments")
cm=Nraw*populationsprint_as_vector(cm, print_symbol=temp_pop_str)
print("\n\t//collision in central moments space")
# print("//calculate equilibrium distributions in cm space")# print("real_t {cm_eq_pop_str}[{q}];\n")# print_as_vector(hardcoded_cm_eq, cm_eq_pop_str) # save time, verbosity# print("//calculate forces in cm space")# print("real_t {F_cm_str}[{q}];")# print_as_vector(hardcoded_F_cm, F_cm_str) # save time, verbosityprint("\t//collide")
defmake_collision(choice):
model_switcher= {
# Relax 2nd moments for hydro, SOI'hydro': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq+ (eye(q) -S_Relax/2) *hardcoded_F_cm,
# Relax 1st moments for ADE, SOI'ade_with_f': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq+ (eye(q) -S_Relax/2) *hardcoded_F_cm,
# Relax 1st moments for ADE, SOI without force'ade': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq,
}
# Get the function from switcher dictionarycm_after_collision=model_switcher.get(choice, lambda: "Invalid argument")
print_as_vector(cm_after_collision, print_symbol=pop_in_str)
make_collision(model)
print("\n\t//back to raw moments")
m=Nraw.inv() *populationsprint_as_vector(m, print_symbol=temp_pop_str)
print("\n\t//back to density-probability functions")
populations=Mraw.inv() *temp_populationsprint_as_vector(populations, print_symbol=pop_in_str)
print("\n}\n")
Output
CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t x_in[9], real_t omega_nu, vector_t u, vector_t F) {
//=== THIS IS AUTOMATICALLY GENERATED CODE ===
real_t uxuy = u.x*u.y;
real_t ux2 = u.x*u.x;
real_t uy2 = u.y*u.y;
real_t m00 = x_in[0] + x_in[1] + x_in[2] + x_in[3] + x_in[4] + x_in[5] + x_in[6] + x_in[7] + x_in[8];
real_t temp[9];for (int i = 0; i < 9; i++) {
temp[i] = x_in[i];}
//raw moments from density-probability functions
x_in[0] = temp[0] + temp[1] + temp[2] + temp[3] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8];
x_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8];
x_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[5] = temp[5] - temp[6] + temp[7] - temp[8];
x_in[6] = temp[5] + temp[6] - temp[7] - temp[8];
x_in[7] = temp[5] - temp[6] - temp[7] + temp[8];
x_in[8] = temp[5] + temp[6] + temp[7] + temp[8];
//central moments from raw moments
temp[0] = x_in[0];
temp[1] = -u.x*x_in[0] + x_in[1];
temp[2] = -u.y*x_in[0] + x_in[2];
temp[3] = ux2*x_in[0] - 2.*u.x*x_in[1] + x_in[3];
temp[4] = uy2*x_in[0] - 2.*u.y*x_in[2] + x_in[4];
temp[5] = uxuy*x_in[0] - u.x*x_in[2] - u.y*x_in[1] + x_in[5];
temp[6] = -ux2*u.y*x_in[0] + ux2*x_in[2] + 2.*uxuy*x_in[1] - 2.*u.x*x_in[5] - u.y*x_in[3] + x_in[6];
temp[7] = -u.x*uy2*x_in[0] + 2.*uxuy*x_in[2] - u.x*x_in[4] + uy2*x_in[1] - 2.*u.y*x_in[5] + x_in[7];
temp[8] = ux2*uy2*x_in[0] - 2.*ux2*u.y*x_in[2] + ux2*x_in[4] - 2.*u.x*uy2*x_in[1] + 4.*uxuy*x_in[5] - 2.*u.x*x_in[7] + uy2*x_in[3] - 2.*u.y*x_in[6] + x_in[8];
//collision in central moments space
//collide
x_in[0] = m00;
x_in[1] = 1/2.*F.x;
x_in[2] = 1/2.*F.y;
x_in[3] = 1/3.*m00*omega_bulk - 1/2.*omega_bulk*temp[3] - 1/2.*omega_bulk*temp[4] - 1/2.*omega_nu*temp[3] + 1/2.*omega_nu*temp[4] + temp[3];
x_in[4] = 1/3.*m00*omega_bulk - 1/2.*omega_bulk*temp[3] - 1/2.*omega_bulk*temp[4] + 1/2.*omega_nu*temp[3] - 1/2.*omega_nu*temp[4] + temp[4];
x_in[5] = -temp[5]*(omega_nu - 1.);
x_in[6] = 1/6.*F.y;
x_in[7] = 1/6.*F.x;
x_in[8] = 1/9.*m00;
//back to raw moments
temp[0] = x_in[0];
temp[1] = u.x*x_in[0] + x_in[1];
temp[2] = u.y*x_in[0] + x_in[2];
temp[3] = ux2*x_in[0] + 2.*u.x*x_in[1] + x_in[3];
temp[4] = uy2*x_in[0] + 2.*u.y*x_in[2] + x_in[4];
temp[5] = uxuy*x_in[0] + u.x*x_in[2] + u.y*x_in[1] + x_in[5];
temp[6] = ux2*u.y*x_in[0] + ux2*x_in[2] + 2.*uxuy*x_in[1] + 2.*u.x*x_in[5] + u.y*x_in[3] + x_in[6];
temp[7] = u.x*uy2*x_in[0] + 2.*uxuy*x_in[2] + u.x*x_in[4] + uy2*x_in[1] + 2.*u.y*x_in[5] + x_in[7];
temp[8] = ux2*uy2*x_in[0] + 2.*ux2*u.y*x_in[2] + ux2*x_in[4] + 2.*u.x*uy2*x_in[1] + 4.*uxuy*x_in[5] + 2.*u.x*x_in[7] + uy2*x_in[3] + 2.*u.y*x_in[6] + x_in[8];
//back to density-probability functions
x_in[0] = temp[0] - temp[3] - temp[4] + temp[8];
x_in[1] = 1/2.*temp[1] + 1/2.*temp[3] - 1/2.*temp[7] - 1/2.*temp[8];
x_in[2] = 1/2.*temp[2] + 1/2.*temp[4] - 1/2.*temp[6] - 1/2.*temp[8];
x_in[3] = -1/2.*temp[1] + 1/2.*temp[3] + 1/2.*temp[7] - 1/2.*temp[8];
x_in[4] = -1/2.*temp[2] + 1/2.*temp[4] + 1/2.*temp[6] - 1/2.*temp[8];
x_in[5] = 1/4.*temp[5] + 1/4.*temp[6] + 1/4.*temp[7] + 1/4.*temp[8];
x_in[6] = -1/4.*temp[5] + 1/4.*temp[6] - 1/4.*temp[7] + 1/4.*temp[8];
x_in[7] = 1/4.*temp[5] - 1/4.*temp[6] - 1/4.*temp[7] + 1/4.*temp[8];
x_in[8] = -1/4.*temp[5] - 1/4.*temp[6] + 1/4.*temp[7] + 1/4.*temp[8];
}
About
Various tools for the TCLB project
Resources
Readme
GPL-3.0 license
Activity
Stars
0 stars
Watchers
1 watching
Forks
0 forks
Report repository
Releases
Packages
Contributors
Languages

, 'i');
if (__m === '*' || __re.test(location.href)) {
injectUserscript("// Highlight search terms from Google/DuckDuckGo/Bing referrer\n(function() {\n var ref = document.referrer;\n var terms = [];\n \n if (ref.includes('google.com') || ref.includes('duckduckgo.com') || ref.includes('bing.com')) {\n var url = new URL(ref);\n var q = url.searchParams.get('q') || url.searchParams.get('p');\n if (q) {\n terms = q.split(/\\s+/).filter(function(t) { return t.length > 2; });\n }\n }\n \n if (terms.length === 0) return;\n \n var style = document.createElement('style');\n style.textContent = '.userscript-highlight { background: #fbbf24; color: #1a1a2e; padding: 1px 3px; border-radius: 2px; }';\n document.head.appendChild(style);\n \n function highlight(node) {\n if (node.nodeType === 3) { // text node\n var text = node.textContent;\n var found = false;\n terms.forEach(function(term) {\n var regex = new RegExp('(' + term.replace(/[.*+?^${}()|[\\]\\\\]/g, '\\\\') + ')', 'gi');\n if (regex.test(text)) {\n found = true;\n var frag = document.createDocumentFragment();\n var parts = text.split(regex);\n parts.forEach(function(part, i) {\n if (i % 2 === 0) {\n frag.appendChild(document.createTextNode(part));\n } else {\n var span = document.createElement('span');\n span.className = 'userscript-highlight';\n span.textContent = part;\n frag.appendChild(span);\n }\n });\n node.parentNode.replaceChild(frag, node);\n }\n });\n } else if (node.nodeType === 1 && node.childNodes) { // element\n var skipTags = ['SCRIPT', 'STYLE', 'NOSCRIPT', 'TEXTAREA', 'INPUT', 'SELECT'];\n if (!skipTags.includes(node.tagName)) {\n Array.from(node.childNodes).forEach(highlight);\n }\n }\n }\n \n highlight(document.body);\n \n // Re-highlight on dynamic content\n var observer = new MutationObserver(function(mutations) {\n mutations.forEach(function(m) {\n m.addedNodes.forEach(function(node) {\n if (node.nodeType === 1 || node.nodeType === 3) highlight(node);\n });\n });\n });\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "Highlight Search Terms");
}
} catch(__e) { console.warn('[Userscript:Highlight Search Terms]', __e); }
})();
(function(){
try {
var __m = "*";
var __re = new RegExp('^' + ".*" + '

Skip to content
BranchesTags
Repository files navigation
TCLB tools
Various tools and extended examples for the TCLB project.
Structure
All tools are divided by language: R, Python, etc. In each directory, the subdirectories should be specific packages related to different aspects of interaction with TCLB
Examples and Papers include scripts used to obtain data used in "real world" cases. Be carefoul - they might be computationally intensive. See Readmi for details
Installation
R
Supervisor: Łukasz Łaniewski-Wołłk
devtools::install_github("CFD-GO/TCLB_tools/R/TCLBtools")
Python
Supervisor: Michał Dzikowski
Content of Python/ directory is meant as a packages repository. Add XXXX/Python to your PYTHONPATH than import packages by subdirectory name.
Python: symbolic_tools
Supervisor: Grzegorz Gruszczyński
It is recomended to use PyCharm to run the scripts.
It facilitates recognition of paths and folders.
Create a PyCharm project by opening the symbolic_tools/ directory in the editor.
Run the examples/ . Some are listed below.
Example: Forcing terms
fromSymbolicCollisions.core.printersimportprint_as_vectorfromSymbolicCollisions.core.ContinousCMTransformsimport \
ContinousCMTransforms, get_mom_vector_from_continuous_deffromSymbolicCollisions.core.cm_symbolsimport \
F3D, dzeta3D, u3DfromSymbolicCollisions.core.cm_symbolsimportrho, moments_dictccmt=ContinousCMTransforms(dzeta3D, u3D, F3D, rho)
lattice='D2Q9'print("\n--- FORCES ---")
print('\n\n// === continuous central moments === \n ')
print('\n//Force -> Force_cm - from continous definition: \n''k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) \n''where fun = forceM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o ')
F_cm=get_mom_vector_from_continuous_def(ccmt.get_force_He_MB,
continuous_transformation=ccmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(F_cm, 'F_cm')
Output:
//Force -> Force_cm - from continous definition: k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) where fun = forceM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o F_cm[0] = 0;
F_cm[1] = Fhydro.x*m00/rho;
F_cm[2] = Fhydro.y*m00/rho;
F_cm[3] = 0;
F_cm[4] = 0;
F_cm[5] = 0;
F_cm[6] = 1/3.*Fhydro.y*m00/rho;
F_cm[7] = 1/3.*Fhydro.x*m00/rho;
F_cm[8] = 0;
Example: Equilibrium Distribution
fromSymbolicCollisions.core.printersimportprint_as_vectorfromsympy.matricesimportMatrixfromsympyimportSymbolfromSymbolicCollisions.core.ContinousCMTransformsimportContinousCMTransforms, get_mom_vector_from_continuous_deffromSymbolicCollisions.core.cm_symbolsimport \
F3D, dzeta3D, u3DfromSymbolicCollisions.core.cm_symbolsimporte_D2Q9, u2D, F2D, rho, moments_dictfromSymbolicCollisions.core.DiscreteCMTransformsimport \
DiscreteCMTransforms, get_mom_vector_from_discrete_deflattice='D2Q9'ccmt=ContinousCMTransforms(dzeta3D, u3D, F3D, rho)
dcmt=DiscreteCMTransforms(e_D2Q9, u2D, F2D, rho)
print("\n--- EQUILIBRIA ---")
print('// === discrete cm ===\n ')
print('\n//population_eq -> cm_eq - by definition: k_mn = sum( (e_ix-ux)^m (e_iy-uy)^n * population_eq_i)')
print("moments: first order (linear) velocity expansion.")
pop_eq=get_mom_vector_from_discrete_def(lambdai: Symbol('m00') *dcmt.get_gamma_first_order(i),
discrete_transform=dcmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(pop_eq, 'cm_eq_first_order')
print("moments: second order (quadratic) velocity expansion.")
pop_eq=get_mom_vector_from_discrete_def(lambdai: Symbol('m00') *dcmt.get_gamma(i),
discrete_transform=dcmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(pop_eq, 'cm_eq_second_order')
print('\n\n// === continous cm === \n ')
# to calculate particular momentrow=moments_dict['D2Q9'][0]
moment=ccmt.get_cm(row, ccmt.get_Maxwellian_DF)
print_as_vector(Matrix([moment]), 'particular_moment')
print('\n//population_eq -> cm_eq - from continous definition: \n''k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) \n''where fun = fM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o ')
cm_eq=get_mom_vector_from_continuous_def(ccmt.get_Maxwellian_DF,
continuous_transformation=ccmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(cm_eq, 'cm_eq')
Output
--- EQUILIBRIA ---
// === discrete cm ===
//population_eq -> cm_eq - by definition: k_mn = sum( (e_ix-ux)^m (e_iy-uy)^n * population_eq_i)
moments: first order (linear) velocity expansion.
cm_eq_first_order[0] = m00;
cm_eq_first_order[1] = 0;
cm_eq_first_order[2] = 0;
cm_eq_first_order[3] = m00*(-ux2 + 1/3.);
cm_eq_first_order[4] = m00*(-uy2 + 1/3.);
cm_eq_first_order[5] = -m00*uxuy;
cm_eq_first_order[6] = 2.*m00*ux2*u.y;
cm_eq_first_order[7] = 2.*m00*u.x*uy2;
cm_eq_first_order[8] = m00*(-3.*ux2*uy2 - 1/3.*ux2 - 1/3.*uy2 + 1/9.);
moments: second order (quadratic) velocity expansion.
cm_eq_second_order[0] = m00;
cm_eq_second_order[1] = 0;
cm_eq_second_order[2] = 0;
cm_eq_second_order[3] = 1/3.*m00;
cm_eq_second_order[4] = 1/3.*m00;
cm_eq_second_order[5] = 0;
cm_eq_second_order[6] = -m00*ux2*u.y;
cm_eq_second_order[7] = -m00*u.x*uy2;
cm_eq_second_order[8] = m00*(3.*ux2*uy2 + 1/9.);
// === continous cm === particular_moment[0] = m00;
//population_eq -> cm_eq - from continous definition: k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) where fun = fM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o cm_eq[0] = m00;
cm_eq[1] = 0;
cm_eq[2] = 0;
cm_eq[3] = 1/3.*m00;
cm_eq[4] = 1/3.*m00;
cm_eq[5] = 0;
cm_eq[6] = 0;
cm_eq[7] = 0;
cm_eq[8] = 1/9.*m00;
Example: Collision Kernel
fromsympy.matricesimporteyefromsympy.printingimportprint_ccodefromSymbolicCollisions.core.cm_symbolsimportomega_ade, omega_b, omega_v, m00fromSymbolicCollisions.core.cm_symbolsimportdynamic_importfromSymbolicCollisions.core.DiscreteCMTransformsimportget_DF, get_m00fromSymbolicCollisions.core.printersimportprint_u2, print_as_vectorfromSymbolicCollisions.core.MatrixGeneratorimportget_raw_moments_matrix, get_shift_matrix# inspired by:# "Consistent Forcing Scheme in the cascaded LBM" L. Fei et al. 2017# eqs 8-12 : (eye(q)-S)*cm + S*cm_eq + (eye(q)-S/2.)*force_in_cm_space# SETUPd=2q=9model='hydro'# choose from '['hydro', 'ade', 'ade_with_f']# DYNAMIC IMPORTSex=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ex_D{d}Q{q}")
ey=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ey_D{d}Q{q}")
ifd==3:
ez=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ez_D{d}Q{q}")
else:
ez=Nonedefget_s_relax_switcher(choice):
s_relax_switcher= {
'hydro': ("SymbolicCollisions.core.cm_symbols", f"S_relax_hydro_D{d}Q{q}"),
'ade_with_f': ("SymbolicCollisions.core.cm_symbols", f"S_relax_ADE_D{d}Q{q}"),
'ade': ("SymbolicCollisions.core.cm_symbols", f"S_relax_ADE_D{d}Q{q}"),
}
which_model=s_relax_switcher.get(choice, lambda: "Invalid argument")
returndynamic_import(*which_model)
S_Relax=get_s_relax_switcher(model)
hardcoded_cm_eq=dynamic_import("SymbolicCollisions.core.hardcoded_results", f"hardcoded_cm_eq_compressible_D{d}Q{q}")
hardcoded_F_cm=dynamic_import("SymbolicCollisions.core.hardcoded_results", f"hardcoded_F_cm_pf_D{d}Q{q}")
fromSymbolicCollisions.core.cm_symbolsimportForce_strasF_str# ARRANGE STUFFMraw=get_raw_moments_matrix(ex, ey, ez)
Nraw=get_shift_matrix(Mraw.inv(), ex, ey, ez)
# from sympy import pprint# pprint(Mraw) # see what you have done# pprint(Nraw)pop_in_str='x_in'# symbol defining populationstemp_pop_str='temp'# symbol defining populationscm_eq_pop_str='cm_eq'# symbol defining populations# GENERATE CODEdefmake_header(choice):
model_switcher= {
'hydro': f"CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t {pop_in_str}[{q}], real_t {omega_v}, vector_t u, vector_t {F_str}) \n{{",
'ade_with_f': f"CudaDeviceFunction void relax_and_collide_ADE_with_F(real_t {pop_in_str}[{q}], real_t {omega_ade}, vector_t u, vector_t {F_str}) \n{{",
'ade': f"CudaDeviceFunction void relax_and_collide_ADE(real_t {pop_in_str}[{q}], real_t {omega_ade}, vector_t u) \n{{",
}
result=model_switcher.get(choice, lambda: "Invalid argument")
print(result)
make_header(model)
print("\t//=== THIS IS AUTOMATICALLY GENERATED CODE ===")
# print(f"real_t {sv} = omega;")# print("real_t bulk_visc = 1./6. ;")# print("real_t {sb} = 1./(3*bulk_visc + 0.5);")# print(f"real_t {sb} = omega_bulk;\n") # s_b = 1./(3*bulk_visc + 0.5)print_u2(d)
print_ccode(get_m00(q, pop_in_str), assign_to=f'\treal_t {m00}')
defmake_variables(choice):
model_switcher= {
'hydro': f"\n\treal_t {temp_pop_str}[{q}];\n",
'ade_with_f': f"\n\treal_t {temp_pop_str}[{q}];\n",
'ade': f"\n\treal_t {temp_pop_str}[{q}];\n",
}
# Get the function from switcher dictionaryresult=model_switcher.get(choice, lambda: "Invalid argument")
print(result)
make_variables(model)
print(f"\tfor (int i = 0; i < {q}; i++) {{\n\t"f"\t{temp_pop_str}[i] = {pop_in_str}[i];}}")
populations=get_DF(q, pop_in_str)
temp_populations=get_DF(q, temp_pop_str)
cm_eq=get_DF(q, cm_eq_pop_str)
F_cm=get_DF(q, F_str)
m=Mraw*temp_populationsprint("\n\t//raw moments from density-probability functions")
# print("\t//[m00, m10, m01, m20, m02, m11, m21, m12, m22]")print_as_vector(m, print_symbol=pop_in_str)
print("\n\t//central moments from raw moments")
cm=Nraw*populationsprint_as_vector(cm, print_symbol=temp_pop_str)
print("\n\t//collision in central moments space")
# print("//calculate equilibrium distributions in cm space")# print("real_t {cm_eq_pop_str}[{q}];\n")# print_as_vector(hardcoded_cm_eq, cm_eq_pop_str) # save time, verbosity# print("//calculate forces in cm space")# print("real_t {F_cm_str}[{q}];")# print_as_vector(hardcoded_F_cm, F_cm_str) # save time, verbosityprint("\t//collide")
defmake_collision(choice):
model_switcher= {
# Relax 2nd moments for hydro, SOI'hydro': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq+ (eye(q) -S_Relax/2) *hardcoded_F_cm,
# Relax 1st moments for ADE, SOI'ade_with_f': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq+ (eye(q) -S_Relax/2) *hardcoded_F_cm,
# Relax 1st moments for ADE, SOI without force'ade': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq,
}
# Get the function from switcher dictionarycm_after_collision=model_switcher.get(choice, lambda: "Invalid argument")
print_as_vector(cm_after_collision, print_symbol=pop_in_str)
make_collision(model)
print("\n\t//back to raw moments")
m=Nraw.inv() *populationsprint_as_vector(m, print_symbol=temp_pop_str)
print("\n\t//back to density-probability functions")
populations=Mraw.inv() *temp_populationsprint_as_vector(populations, print_symbol=pop_in_str)
print("\n}\n")
Output
CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t x_in[9], real_t omega_nu, vector_t u, vector_t F) {
//=== THIS IS AUTOMATICALLY GENERATED CODE ===
real_t uxuy = u.x*u.y;
real_t ux2 = u.x*u.x;
real_t uy2 = u.y*u.y;
real_t m00 = x_in[0] + x_in[1] + x_in[2] + x_in[3] + x_in[4] + x_in[5] + x_in[6] + x_in[7] + x_in[8];
real_t temp[9];for (int i = 0; i < 9; i++) {
temp[i] = x_in[i];}
//raw moments from density-probability functions
x_in[0] = temp[0] + temp[1] + temp[2] + temp[3] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8];
x_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8];
x_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[5] = temp[5] - temp[6] + temp[7] - temp[8];
x_in[6] = temp[5] + temp[6] - temp[7] - temp[8];
x_in[7] = temp[5] - temp[6] - temp[7] + temp[8];
x_in[8] = temp[5] + temp[6] + temp[7] + temp[8];
//central moments from raw moments
temp[0] = x_in[0];
temp[1] = -u.x*x_in[0] + x_in[1];
temp[2] = -u.y*x_in[0] + x_in[2];
temp[3] = ux2*x_in[0] - 2.*u.x*x_in[1] + x_in[3];
temp[4] = uy2*x_in[0] - 2.*u.y*x_in[2] + x_in[4];
temp[5] = uxuy*x_in[0] - u.x*x_in[2] - u.y*x_in[1] + x_in[5];
temp[6] = -ux2*u.y*x_in[0] + ux2*x_in[2] + 2.*uxuy*x_in[1] - 2.*u.x*x_in[5] - u.y*x_in[3] + x_in[6];
temp[7] = -u.x*uy2*x_in[0] + 2.*uxuy*x_in[2] - u.x*x_in[4] + uy2*x_in[1] - 2.*u.y*x_in[5] + x_in[7];
temp[8] = ux2*uy2*x_in[0] - 2.*ux2*u.y*x_in[2] + ux2*x_in[4] - 2.*u.x*uy2*x_in[1] + 4.*uxuy*x_in[5] - 2.*u.x*x_in[7] + uy2*x_in[3] - 2.*u.y*x_in[6] + x_in[8];
//collision in central moments space
//collide
x_in[0] = m00;
x_in[1] = 1/2.*F.x;
x_in[2] = 1/2.*F.y;
x_in[3] = 1/3.*m00*omega_bulk - 1/2.*omega_bulk*temp[3] - 1/2.*omega_bulk*temp[4] - 1/2.*omega_nu*temp[3] + 1/2.*omega_nu*temp[4] + temp[3];
x_in[4] = 1/3.*m00*omega_bulk - 1/2.*omega_bulk*temp[3] - 1/2.*omega_bulk*temp[4] + 1/2.*omega_nu*temp[3] - 1/2.*omega_nu*temp[4] + temp[4];
x_in[5] = -temp[5]*(omega_nu - 1.);
x_in[6] = 1/6.*F.y;
x_in[7] = 1/6.*F.x;
x_in[8] = 1/9.*m00;
//back to raw moments
temp[0] = x_in[0];
temp[1] = u.x*x_in[0] + x_in[1];
temp[2] = u.y*x_in[0] + x_in[2];
temp[3] = ux2*x_in[0] + 2.*u.x*x_in[1] + x_in[3];
temp[4] = uy2*x_in[0] + 2.*u.y*x_in[2] + x_in[4];
temp[5] = uxuy*x_in[0] + u.x*x_in[2] + u.y*x_in[1] + x_in[5];
temp[6] = ux2*u.y*x_in[0] + ux2*x_in[2] + 2.*uxuy*x_in[1] + 2.*u.x*x_in[5] + u.y*x_in[3] + x_in[6];
temp[7] = u.x*uy2*x_in[0] + 2.*uxuy*x_in[2] + u.x*x_in[4] + uy2*x_in[1] + 2.*u.y*x_in[5] + x_in[7];
temp[8] = ux2*uy2*x_in[0] + 2.*ux2*u.y*x_in[2] + ux2*x_in[4] + 2.*u.x*uy2*x_in[1] + 4.*uxuy*x_in[5] + 2.*u.x*x_in[7] + uy2*x_in[3] + 2.*u.y*x_in[6] + x_in[8];
//back to density-probability functions
x_in[0] = temp[0] - temp[3] - temp[4] + temp[8];
x_in[1] = 1/2.*temp[1] + 1/2.*temp[3] - 1/2.*temp[7] - 1/2.*temp[8];
x_in[2] = 1/2.*temp[2] + 1/2.*temp[4] - 1/2.*temp[6] - 1/2.*temp[8];
x_in[3] = -1/2.*temp[1] + 1/2.*temp[3] + 1/2.*temp[7] - 1/2.*temp[8];
x_in[4] = -1/2.*temp[2] + 1/2.*temp[4] + 1/2.*temp[6] - 1/2.*temp[8];
x_in[5] = 1/4.*temp[5] + 1/4.*temp[6] + 1/4.*temp[7] + 1/4.*temp[8];
x_in[6] = -1/4.*temp[5] + 1/4.*temp[6] - 1/4.*temp[7] + 1/4.*temp[8];
x_in[7] = 1/4.*temp[5] - 1/4.*temp[6] - 1/4.*temp[7] + 1/4.*temp[8];
x_in[8] = -1/4.*temp[5] - 1/4.*temp[6] + 1/4.*temp[7] + 1/4.*temp[8];
}
About
Various tools for the TCLB project
Resources
Readme
GPL-3.0 license
Activity
Stars
0 stars
Watchers
1 watching
Forks
0 forks
Report repository
Releases
Packages
Contributors
Languages

, 'i');
if (__m === '*' || __re.test(location.href)) {
injectUserscript("// Strip utm_, fbclid, gclid, etc. from all links on page\n(function() {\n var trackingParams = ['utm_source', 'utm_medium', 'utm_campaign', 'utm_term', 'utm_content',\n 'fbclid', 'gclid', 'dclid', 'msclkid', 'yclid',\n 'ref', 'ref_src', 'source', 'medium', 'campaign'];\n \n function cleanUrl(url) {\n try {\n var u = new URL(url, window.location.origin);\n var changed = false;\n trackingParams.forEach(function(p) {\n if (u.searchParams.has(p)) {\n u.searchParams.delete(p);\n changed = true;\n }\n });\n return changed ? u.toString() : url;\n } catch (e) {\n return url;\n }\n }\n \n function cleanLinks() {\n document.querySelectorAll('a[href]').forEach(function(a) {\n var clean = cleanUrl(a.href);\n if (clean !== a.href) a.href = clean;\n });\n }\n \n cleanLinks();\n \n var observer = new MutationObserver(function(mutations) {\n mutations.forEach(function(m) {\n m.addedNodes.forEach(function(node) {\n if (node.nodeType === 1) {\n if (node.tagName === 'A') cleanLinks();\n node.querySelectorAll('a[href]').forEach(function(a) {\n var clean = cleanUrl(a.href);\n if (clean !== a.href) a.href = clean;\n });\n }\n });\n });\n });\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "Remove Tracking Parameters from Links");
}
} catch(__e) { console.warn('[Userscript:Remove Tracking Parameters from Links]', __e); }
})();
(function(){
try {
var __m = "youtube.com";
var __re = new RegExp('^' + "youtube\\.com" + '

Skip to content
BranchesTags
Repository files navigation
TCLB tools
Various tools and extended examples for the TCLB project.
Structure
All tools are divided by language: R, Python, etc. In each directory, the subdirectories should be specific packages related to different aspects of interaction with TCLB
Examples and Papers include scripts used to obtain data used in "real world" cases. Be carefoul - they might be computationally intensive. See Readmi for details
Installation
R
Supervisor: Łukasz Łaniewski-Wołłk
devtools::install_github("CFD-GO/TCLB_tools/R/TCLBtools")
Python
Supervisor: Michał Dzikowski
Content of Python/ directory is meant as a packages repository. Add XXXX/Python to your PYTHONPATH than import packages by subdirectory name.
Python: symbolic_tools
Supervisor: Grzegorz Gruszczyński
It is recomended to use PyCharm to run the scripts.
It facilitates recognition of paths and folders.
Create a PyCharm project by opening the symbolic_tools/ directory in the editor.
Run the examples/ . Some are listed below.
Example: Forcing terms
fromSymbolicCollisions.core.printersimportprint_as_vectorfromSymbolicCollisions.core.ContinousCMTransformsimport \
ContinousCMTransforms, get_mom_vector_from_continuous_deffromSymbolicCollisions.core.cm_symbolsimport \
F3D, dzeta3D, u3DfromSymbolicCollisions.core.cm_symbolsimportrho, moments_dictccmt=ContinousCMTransforms(dzeta3D, u3D, F3D, rho)
lattice='D2Q9'print("\n--- FORCES ---")
print('\n\n// === continuous central moments === \n ')
print('\n//Force -> Force_cm - from continous definition: \n''k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) \n''where fun = forceM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o ')
F_cm=get_mom_vector_from_continuous_def(ccmt.get_force_He_MB,
continuous_transformation=ccmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(F_cm, 'F_cm')
Output:
//Force -> Force_cm - from continous definition: k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) where fun = forceM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o F_cm[0] = 0;
F_cm[1] = Fhydro.x*m00/rho;
F_cm[2] = Fhydro.y*m00/rho;
F_cm[3] = 0;
F_cm[4] = 0;
F_cm[5] = 0;
F_cm[6] = 1/3.*Fhydro.y*m00/rho;
F_cm[7] = 1/3.*Fhydro.x*m00/rho;
F_cm[8] = 0;
Example: Equilibrium Distribution
fromSymbolicCollisions.core.printersimportprint_as_vectorfromsympy.matricesimportMatrixfromsympyimportSymbolfromSymbolicCollisions.core.ContinousCMTransformsimportContinousCMTransforms, get_mom_vector_from_continuous_deffromSymbolicCollisions.core.cm_symbolsimport \
F3D, dzeta3D, u3DfromSymbolicCollisions.core.cm_symbolsimporte_D2Q9, u2D, F2D, rho, moments_dictfromSymbolicCollisions.core.DiscreteCMTransformsimport \
DiscreteCMTransforms, get_mom_vector_from_discrete_deflattice='D2Q9'ccmt=ContinousCMTransforms(dzeta3D, u3D, F3D, rho)
dcmt=DiscreteCMTransforms(e_D2Q9, u2D, F2D, rho)
print("\n--- EQUILIBRIA ---")
print('// === discrete cm ===\n ')
print('\n//population_eq -> cm_eq - by definition: k_mn = sum( (e_ix-ux)^m (e_iy-uy)^n * population_eq_i)')
print("moments: first order (linear) velocity expansion.")
pop_eq=get_mom_vector_from_discrete_def(lambdai: Symbol('m00') *dcmt.get_gamma_first_order(i),
discrete_transform=dcmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(pop_eq, 'cm_eq_first_order')
print("moments: second order (quadratic) velocity expansion.")
pop_eq=get_mom_vector_from_discrete_def(lambdai: Symbol('m00') *dcmt.get_gamma(i),
discrete_transform=dcmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(pop_eq, 'cm_eq_second_order')
print('\n\n// === continous cm === \n ')
# to calculate particular momentrow=moments_dict['D2Q9'][0]
moment=ccmt.get_cm(row, ccmt.get_Maxwellian_DF)
print_as_vector(Matrix([moment]), 'particular_moment')
print('\n//population_eq -> cm_eq - from continous definition: \n''k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) \n''where fun = fM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o ')
cm_eq=get_mom_vector_from_continuous_def(ccmt.get_Maxwellian_DF,
continuous_transformation=ccmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(cm_eq, 'cm_eq')
Output
--- EQUILIBRIA ---
// === discrete cm ===
//population_eq -> cm_eq - by definition: k_mn = sum( (e_ix-ux)^m (e_iy-uy)^n * population_eq_i)
moments: first order (linear) velocity expansion.
cm_eq_first_order[0] = m00;
cm_eq_first_order[1] = 0;
cm_eq_first_order[2] = 0;
cm_eq_first_order[3] = m00*(-ux2 + 1/3.);
cm_eq_first_order[4] = m00*(-uy2 + 1/3.);
cm_eq_first_order[5] = -m00*uxuy;
cm_eq_first_order[6] = 2.*m00*ux2*u.y;
cm_eq_first_order[7] = 2.*m00*u.x*uy2;
cm_eq_first_order[8] = m00*(-3.*ux2*uy2 - 1/3.*ux2 - 1/3.*uy2 + 1/9.);
moments: second order (quadratic) velocity expansion.
cm_eq_second_order[0] = m00;
cm_eq_second_order[1] = 0;
cm_eq_second_order[2] = 0;
cm_eq_second_order[3] = 1/3.*m00;
cm_eq_second_order[4] = 1/3.*m00;
cm_eq_second_order[5] = 0;
cm_eq_second_order[6] = -m00*ux2*u.y;
cm_eq_second_order[7] = -m00*u.x*uy2;
cm_eq_second_order[8] = m00*(3.*ux2*uy2 + 1/9.);
// === continous cm === particular_moment[0] = m00;
//population_eq -> cm_eq - from continous definition: k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) where fun = fM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o cm_eq[0] = m00;
cm_eq[1] = 0;
cm_eq[2] = 0;
cm_eq[3] = 1/3.*m00;
cm_eq[4] = 1/3.*m00;
cm_eq[5] = 0;
cm_eq[6] = 0;
cm_eq[7] = 0;
cm_eq[8] = 1/9.*m00;
Example: Collision Kernel
fromsympy.matricesimporteyefromsympy.printingimportprint_ccodefromSymbolicCollisions.core.cm_symbolsimportomega_ade, omega_b, omega_v, m00fromSymbolicCollisions.core.cm_symbolsimportdynamic_importfromSymbolicCollisions.core.DiscreteCMTransformsimportget_DF, get_m00fromSymbolicCollisions.core.printersimportprint_u2, print_as_vectorfromSymbolicCollisions.core.MatrixGeneratorimportget_raw_moments_matrix, get_shift_matrix# inspired by:# "Consistent Forcing Scheme in the cascaded LBM" L. Fei et al. 2017# eqs 8-12 : (eye(q)-S)*cm + S*cm_eq + (eye(q)-S/2.)*force_in_cm_space# SETUPd=2q=9model='hydro'# choose from '['hydro', 'ade', 'ade_with_f']# DYNAMIC IMPORTSex=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ex_D{d}Q{q}")
ey=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ey_D{d}Q{q}")
ifd==3:
ez=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ez_D{d}Q{q}")
else:
ez=Nonedefget_s_relax_switcher(choice):
s_relax_switcher= {
'hydro': ("SymbolicCollisions.core.cm_symbols", f"S_relax_hydro_D{d}Q{q}"),
'ade_with_f': ("SymbolicCollisions.core.cm_symbols", f"S_relax_ADE_D{d}Q{q}"),
'ade': ("SymbolicCollisions.core.cm_symbols", f"S_relax_ADE_D{d}Q{q}"),
}
which_model=s_relax_switcher.get(choice, lambda: "Invalid argument")
returndynamic_import(*which_model)
S_Relax=get_s_relax_switcher(model)
hardcoded_cm_eq=dynamic_import("SymbolicCollisions.core.hardcoded_results", f"hardcoded_cm_eq_compressible_D{d}Q{q}")
hardcoded_F_cm=dynamic_import("SymbolicCollisions.core.hardcoded_results", f"hardcoded_F_cm_pf_D{d}Q{q}")
fromSymbolicCollisions.core.cm_symbolsimportForce_strasF_str# ARRANGE STUFFMraw=get_raw_moments_matrix(ex, ey, ez)
Nraw=get_shift_matrix(Mraw.inv(), ex, ey, ez)
# from sympy import pprint# pprint(Mraw) # see what you have done# pprint(Nraw)pop_in_str='x_in'# symbol defining populationstemp_pop_str='temp'# symbol defining populationscm_eq_pop_str='cm_eq'# symbol defining populations# GENERATE CODEdefmake_header(choice):
model_switcher= {
'hydro': f"CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t {pop_in_str}[{q}], real_t {omega_v}, vector_t u, vector_t {F_str}) \n{{",
'ade_with_f': f"CudaDeviceFunction void relax_and_collide_ADE_with_F(real_t {pop_in_str}[{q}], real_t {omega_ade}, vector_t u, vector_t {F_str}) \n{{",
'ade': f"CudaDeviceFunction void relax_and_collide_ADE(real_t {pop_in_str}[{q}], real_t {omega_ade}, vector_t u) \n{{",
}
result=model_switcher.get(choice, lambda: "Invalid argument")
print(result)
make_header(model)
print("\t//=== THIS IS AUTOMATICALLY GENERATED CODE ===")
# print(f"real_t {sv} = omega;")# print("real_t bulk_visc = 1./6. ;")# print("real_t {sb} = 1./(3*bulk_visc + 0.5);")# print(f"real_t {sb} = omega_bulk;\n") # s_b = 1./(3*bulk_visc + 0.5)print_u2(d)
print_ccode(get_m00(q, pop_in_str), assign_to=f'\treal_t {m00}')
defmake_variables(choice):
model_switcher= {
'hydro': f"\n\treal_t {temp_pop_str}[{q}];\n",
'ade_with_f': f"\n\treal_t {temp_pop_str}[{q}];\n",
'ade': f"\n\treal_t {temp_pop_str}[{q}];\n",
}
# Get the function from switcher dictionaryresult=model_switcher.get(choice, lambda: "Invalid argument")
print(result)
make_variables(model)
print(f"\tfor (int i = 0; i < {q}; i++) {{\n\t"f"\t{temp_pop_str}[i] = {pop_in_str}[i];}}")
populations=get_DF(q, pop_in_str)
temp_populations=get_DF(q, temp_pop_str)
cm_eq=get_DF(q, cm_eq_pop_str)
F_cm=get_DF(q, F_str)
m=Mraw*temp_populationsprint("\n\t//raw moments from density-probability functions")
# print("\t//[m00, m10, m01, m20, m02, m11, m21, m12, m22]")print_as_vector(m, print_symbol=pop_in_str)
print("\n\t//central moments from raw moments")
cm=Nraw*populationsprint_as_vector(cm, print_symbol=temp_pop_str)
print("\n\t//collision in central moments space")
# print("//calculate equilibrium distributions in cm space")# print("real_t {cm_eq_pop_str}[{q}];\n")# print_as_vector(hardcoded_cm_eq, cm_eq_pop_str) # save time, verbosity# print("//calculate forces in cm space")# print("real_t {F_cm_str}[{q}];")# print_as_vector(hardcoded_F_cm, F_cm_str) # save time, verbosityprint("\t//collide")
defmake_collision(choice):
model_switcher= {
# Relax 2nd moments for hydro, SOI'hydro': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq+ (eye(q) -S_Relax/2) *hardcoded_F_cm,
# Relax 1st moments for ADE, SOI'ade_with_f': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq+ (eye(q) -S_Relax/2) *hardcoded_F_cm,
# Relax 1st moments for ADE, SOI without force'ade': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq,
}
# Get the function from switcher dictionarycm_after_collision=model_switcher.get(choice, lambda: "Invalid argument")
print_as_vector(cm_after_collision, print_symbol=pop_in_str)
make_collision(model)
print("\n\t//back to raw moments")
m=Nraw.inv() *populationsprint_as_vector(m, print_symbol=temp_pop_str)
print("\n\t//back to density-probability functions")
populations=Mraw.inv() *temp_populationsprint_as_vector(populations, print_symbol=pop_in_str)
print("\n}\n")
Output
CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t x_in[9], real_t omega_nu, vector_t u, vector_t F) {
//=== THIS IS AUTOMATICALLY GENERATED CODE ===
real_t uxuy = u.x*u.y;
real_t ux2 = u.x*u.x;
real_t uy2 = u.y*u.y;
real_t m00 = x_in[0] + x_in[1] + x_in[2] + x_in[3] + x_in[4] + x_in[5] + x_in[6] + x_in[7] + x_in[8];
real_t temp[9];for (int i = 0; i < 9; i++) {
temp[i] = x_in[i];}
//raw moments from density-probability functions
x_in[0] = temp[0] + temp[1] + temp[2] + temp[3] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8];
x_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8];
x_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[5] = temp[5] - temp[6] + temp[7] - temp[8];
x_in[6] = temp[5] + temp[6] - temp[7] - temp[8];
x_in[7] = temp[5] - temp[6] - temp[7] + temp[8];
x_in[8] = temp[5] + temp[6] + temp[7] + temp[8];
//central moments from raw moments
temp[0] = x_in[0];
temp[1] = -u.x*x_in[0] + x_in[1];
temp[2] = -u.y*x_in[0] + x_in[2];
temp[3] = ux2*x_in[0] - 2.*u.x*x_in[1] + x_in[3];
temp[4] = uy2*x_in[0] - 2.*u.y*x_in[2] + x_in[4];
temp[5] = uxuy*x_in[0] - u.x*x_in[2] - u.y*x_in[1] + x_in[5];
temp[6] = -ux2*u.y*x_in[0] + ux2*x_in[2] + 2.*uxuy*x_in[1] - 2.*u.x*x_in[5] - u.y*x_in[3] + x_in[6];
temp[7] = -u.x*uy2*x_in[0] + 2.*uxuy*x_in[2] - u.x*x_in[4] + uy2*x_in[1] - 2.*u.y*x_in[5] + x_in[7];
temp[8] = ux2*uy2*x_in[0] - 2.*ux2*u.y*x_in[2] + ux2*x_in[4] - 2.*u.x*uy2*x_in[1] + 4.*uxuy*x_in[5] - 2.*u.x*x_in[7] + uy2*x_in[3] - 2.*u.y*x_in[6] + x_in[8];
//collision in central moments space
//collide
x_in[0] = m00;
x_in[1] = 1/2.*F.x;
x_in[2] = 1/2.*F.y;
x_in[3] = 1/3.*m00*omega_bulk - 1/2.*omega_bulk*temp[3] - 1/2.*omega_bulk*temp[4] - 1/2.*omega_nu*temp[3] + 1/2.*omega_nu*temp[4] + temp[3];
x_in[4] = 1/3.*m00*omega_bulk - 1/2.*omega_bulk*temp[3] - 1/2.*omega_bulk*temp[4] + 1/2.*omega_nu*temp[3] - 1/2.*omega_nu*temp[4] + temp[4];
x_in[5] = -temp[5]*(omega_nu - 1.);
x_in[6] = 1/6.*F.y;
x_in[7] = 1/6.*F.x;
x_in[8] = 1/9.*m00;
//back to raw moments
temp[0] = x_in[0];
temp[1] = u.x*x_in[0] + x_in[1];
temp[2] = u.y*x_in[0] + x_in[2];
temp[3] = ux2*x_in[0] + 2.*u.x*x_in[1] + x_in[3];
temp[4] = uy2*x_in[0] + 2.*u.y*x_in[2] + x_in[4];
temp[5] = uxuy*x_in[0] + u.x*x_in[2] + u.y*x_in[1] + x_in[5];
temp[6] = ux2*u.y*x_in[0] + ux2*x_in[2] + 2.*uxuy*x_in[1] + 2.*u.x*x_in[5] + u.y*x_in[3] + x_in[6];
temp[7] = u.x*uy2*x_in[0] + 2.*uxuy*x_in[2] + u.x*x_in[4] + uy2*x_in[1] + 2.*u.y*x_in[5] + x_in[7];
temp[8] = ux2*uy2*x_in[0] + 2.*ux2*u.y*x_in[2] + ux2*x_in[4] + 2.*u.x*uy2*x_in[1] + 4.*uxuy*x_in[5] + 2.*u.x*x_in[7] + uy2*x_in[3] + 2.*u.y*x_in[6] + x_in[8];
//back to density-probability functions
x_in[0] = temp[0] - temp[3] - temp[4] + temp[8];
x_in[1] = 1/2.*temp[1] + 1/2.*temp[3] - 1/2.*temp[7] - 1/2.*temp[8];
x_in[2] = 1/2.*temp[2] + 1/2.*temp[4] - 1/2.*temp[6] - 1/2.*temp[8];
x_in[3] = -1/2.*temp[1] + 1/2.*temp[3] + 1/2.*temp[7] - 1/2.*temp[8];
x_in[4] = -1/2.*temp[2] + 1/2.*temp[4] + 1/2.*temp[6] - 1/2.*temp[8];
x_in[5] = 1/4.*temp[5] + 1/4.*temp[6] + 1/4.*temp[7] + 1/4.*temp[8];
x_in[6] = -1/4.*temp[5] + 1/4.*temp[6] - 1/4.*temp[7] + 1/4.*temp[8];
x_in[7] = 1/4.*temp[5] - 1/4.*temp[6] - 1/4.*temp[7] + 1/4.*temp[8];
x_in[8] = -1/4.*temp[5] - 1/4.*temp[6] + 1/4.*temp[7] + 1/4.*temp[8];
}
About
Various tools for the TCLB project
Resources
Readme
GPL-3.0 license
Activity
Stars
0 stars
Watchers
1 watching
Forks
0 forks
Report repository
Releases
Packages
Contributors
Languages

, 'i');
if (__m === '*' || __re.test(location.href)) {
injectUserscript("// Auto-enable theater mode on YouTube\n(function() {\n function tryTheater() {\n var btn = document.querySelector('button[aria-label=\"Theater mode\"], ytd-player #player button[title=\"Theater mode\"]');\n if (btn && !btn.classList.contains('activated')) {\n btn.click();\n }\n }\n \n // Try immediately\n tryTheater();\n \n // Try after navigation (SPA)\n var lastUrl = location.href;\n setInterval(function() {\n if (location.href !== lastUrl) {\n lastUrl = location.href;\n setTimeout(tryTheater, 500);\n }\n }, 1000);\n \n // Also try on player load\n var observer = new MutationObserver(tryTheater);\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "YouTube Theater Mode Default");
}
} catch(__e) { console.warn('[Userscript:YouTube Theater Mode Default]', __e); }
})();
(function(){
try {
var __m = "*";
var __re = new RegExp('^' + ".*" + '

Skip to content
BranchesTags
Repository files navigation
TCLB tools
Various tools and extended examples for the TCLB project.
Structure
All tools are divided by language: R, Python, etc. In each directory, the subdirectories should be specific packages related to different aspects of interaction with TCLB
Examples and Papers include scripts used to obtain data used in "real world" cases. Be carefoul - they might be computationally intensive. See Readmi for details
Installation
R
Supervisor: Łukasz Łaniewski-Wołłk
devtools::install_github("CFD-GO/TCLB_tools/R/TCLBtools")
Python
Supervisor: Michał Dzikowski
Content of Python/ directory is meant as a packages repository. Add XXXX/Python to your PYTHONPATH than import packages by subdirectory name.
Python: symbolic_tools
Supervisor: Grzegorz Gruszczyński
It is recomended to use PyCharm to run the scripts.
It facilitates recognition of paths and folders.
Create a PyCharm project by opening the symbolic_tools/ directory in the editor.
Run the examples/ . Some are listed below.
Example: Forcing terms
fromSymbolicCollisions.core.printersimportprint_as_vectorfromSymbolicCollisions.core.ContinousCMTransformsimport \
ContinousCMTransforms, get_mom_vector_from_continuous_deffromSymbolicCollisions.core.cm_symbolsimport \
F3D, dzeta3D, u3DfromSymbolicCollisions.core.cm_symbolsimportrho, moments_dictccmt=ContinousCMTransforms(dzeta3D, u3D, F3D, rho)
lattice='D2Q9'print("\n--- FORCES ---")
print('\n\n// === continuous central moments === \n ')
print('\n//Force -> Force_cm - from continous definition: \n''k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) \n''where fun = forceM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o ')
F_cm=get_mom_vector_from_continuous_def(ccmt.get_force_He_MB,
continuous_transformation=ccmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(F_cm, 'F_cm')
Output:
//Force -> Force_cm - from continous definition: k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) where fun = forceM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o F_cm[0] = 0;
F_cm[1] = Fhydro.x*m00/rho;
F_cm[2] = Fhydro.y*m00/rho;
F_cm[3] = 0;
F_cm[4] = 0;
F_cm[5] = 0;
F_cm[6] = 1/3.*Fhydro.y*m00/rho;
F_cm[7] = 1/3.*Fhydro.x*m00/rho;
F_cm[8] = 0;
Example: Equilibrium Distribution
fromSymbolicCollisions.core.printersimportprint_as_vectorfromsympy.matricesimportMatrixfromsympyimportSymbolfromSymbolicCollisions.core.ContinousCMTransformsimportContinousCMTransforms, get_mom_vector_from_continuous_deffromSymbolicCollisions.core.cm_symbolsimport \
F3D, dzeta3D, u3DfromSymbolicCollisions.core.cm_symbolsimporte_D2Q9, u2D, F2D, rho, moments_dictfromSymbolicCollisions.core.DiscreteCMTransformsimport \
DiscreteCMTransforms, get_mom_vector_from_discrete_deflattice='D2Q9'ccmt=ContinousCMTransforms(dzeta3D, u3D, F3D, rho)
dcmt=DiscreteCMTransforms(e_D2Q9, u2D, F2D, rho)
print("\n--- EQUILIBRIA ---")
print('// === discrete cm ===\n ')
print('\n//population_eq -> cm_eq - by definition: k_mn = sum( (e_ix-ux)^m (e_iy-uy)^n * population_eq_i)')
print("moments: first order (linear) velocity expansion.")
pop_eq=get_mom_vector_from_discrete_def(lambdai: Symbol('m00') *dcmt.get_gamma_first_order(i),
discrete_transform=dcmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(pop_eq, 'cm_eq_first_order')
print("moments: second order (quadratic) velocity expansion.")
pop_eq=get_mom_vector_from_discrete_def(lambdai: Symbol('m00') *dcmt.get_gamma(i),
discrete_transform=dcmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(pop_eq, 'cm_eq_second_order')
print('\n\n// === continous cm === \n ')
# to calculate particular momentrow=moments_dict['D2Q9'][0]
moment=ccmt.get_cm(row, ccmt.get_Maxwellian_DF)
print_as_vector(Matrix([moment]), 'particular_moment')
print('\n//population_eq -> cm_eq - from continous definition: \n''k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) \n''where fun = fM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o ')
cm_eq=get_mom_vector_from_continuous_def(ccmt.get_Maxwellian_DF,
continuous_transformation=ccmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(cm_eq, 'cm_eq')
Output
--- EQUILIBRIA ---
// === discrete cm ===
//population_eq -> cm_eq - by definition: k_mn = sum( (e_ix-ux)^m (e_iy-uy)^n * population_eq_i)
moments: first order (linear) velocity expansion.
cm_eq_first_order[0] = m00;
cm_eq_first_order[1] = 0;
cm_eq_first_order[2] = 0;
cm_eq_first_order[3] = m00*(-ux2 + 1/3.);
cm_eq_first_order[4] = m00*(-uy2 + 1/3.);
cm_eq_first_order[5] = -m00*uxuy;
cm_eq_first_order[6] = 2.*m00*ux2*u.y;
cm_eq_first_order[7] = 2.*m00*u.x*uy2;
cm_eq_first_order[8] = m00*(-3.*ux2*uy2 - 1/3.*ux2 - 1/3.*uy2 + 1/9.);
moments: second order (quadratic) velocity expansion.
cm_eq_second_order[0] = m00;
cm_eq_second_order[1] = 0;
cm_eq_second_order[2] = 0;
cm_eq_second_order[3] = 1/3.*m00;
cm_eq_second_order[4] = 1/3.*m00;
cm_eq_second_order[5] = 0;
cm_eq_second_order[6] = -m00*ux2*u.y;
cm_eq_second_order[7] = -m00*u.x*uy2;
cm_eq_second_order[8] = m00*(3.*ux2*uy2 + 1/9.);
// === continous cm === particular_moment[0] = m00;
//population_eq -> cm_eq - from continous definition: k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) where fun = fM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o cm_eq[0] = m00;
cm_eq[1] = 0;
cm_eq[2] = 0;
cm_eq[3] = 1/3.*m00;
cm_eq[4] = 1/3.*m00;
cm_eq[5] = 0;
cm_eq[6] = 0;
cm_eq[7] = 0;
cm_eq[8] = 1/9.*m00;
Example: Collision Kernel
fromsympy.matricesimporteyefromsympy.printingimportprint_ccodefromSymbolicCollisions.core.cm_symbolsimportomega_ade, omega_b, omega_v, m00fromSymbolicCollisions.core.cm_symbolsimportdynamic_importfromSymbolicCollisions.core.DiscreteCMTransformsimportget_DF, get_m00fromSymbolicCollisions.core.printersimportprint_u2, print_as_vectorfromSymbolicCollisions.core.MatrixGeneratorimportget_raw_moments_matrix, get_shift_matrix# inspired by:# "Consistent Forcing Scheme in the cascaded LBM" L. Fei et al. 2017# eqs 8-12 : (eye(q)-S)*cm + S*cm_eq + (eye(q)-S/2.)*force_in_cm_space# SETUPd=2q=9model='hydro'# choose from '['hydro', 'ade', 'ade_with_f']# DYNAMIC IMPORTSex=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ex_D{d}Q{q}")
ey=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ey_D{d}Q{q}")
ifd==3:
ez=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ez_D{d}Q{q}")
else:
ez=Nonedefget_s_relax_switcher(choice):
s_relax_switcher= {
'hydro': ("SymbolicCollisions.core.cm_symbols", f"S_relax_hydro_D{d}Q{q}"),
'ade_with_f': ("SymbolicCollisions.core.cm_symbols", f"S_relax_ADE_D{d}Q{q}"),
'ade': ("SymbolicCollisions.core.cm_symbols", f"S_relax_ADE_D{d}Q{q}"),
}
which_model=s_relax_switcher.get(choice, lambda: "Invalid argument")
returndynamic_import(*which_model)
S_Relax=get_s_relax_switcher(model)
hardcoded_cm_eq=dynamic_import("SymbolicCollisions.core.hardcoded_results", f"hardcoded_cm_eq_compressible_D{d}Q{q}")
hardcoded_F_cm=dynamic_import("SymbolicCollisions.core.hardcoded_results", f"hardcoded_F_cm_pf_D{d}Q{q}")
fromSymbolicCollisions.core.cm_symbolsimportForce_strasF_str# ARRANGE STUFFMraw=get_raw_moments_matrix(ex, ey, ez)
Nraw=get_shift_matrix(Mraw.inv(), ex, ey, ez)
# from sympy import pprint# pprint(Mraw) # see what you have done# pprint(Nraw)pop_in_str='x_in'# symbol defining populationstemp_pop_str='temp'# symbol defining populationscm_eq_pop_str='cm_eq'# symbol defining populations# GENERATE CODEdefmake_header(choice):
model_switcher= {
'hydro': f"CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t {pop_in_str}[{q}], real_t {omega_v}, vector_t u, vector_t {F_str}) \n{{",
'ade_with_f': f"CudaDeviceFunction void relax_and_collide_ADE_with_F(real_t {pop_in_str}[{q}], real_t {omega_ade}, vector_t u, vector_t {F_str}) \n{{",
'ade': f"CudaDeviceFunction void relax_and_collide_ADE(real_t {pop_in_str}[{q}], real_t {omega_ade}, vector_t u) \n{{",
}
result=model_switcher.get(choice, lambda: "Invalid argument")
print(result)
make_header(model)
print("\t//=== THIS IS AUTOMATICALLY GENERATED CODE ===")
# print(f"real_t {sv} = omega;")# print("real_t bulk_visc = 1./6. ;")# print("real_t {sb} = 1./(3*bulk_visc + 0.5);")# print(f"real_t {sb} = omega_bulk;\n") # s_b = 1./(3*bulk_visc + 0.5)print_u2(d)
print_ccode(get_m00(q, pop_in_str), assign_to=f'\treal_t {m00}')
defmake_variables(choice):
model_switcher= {
'hydro': f"\n\treal_t {temp_pop_str}[{q}];\n",
'ade_with_f': f"\n\treal_t {temp_pop_str}[{q}];\n",
'ade': f"\n\treal_t {temp_pop_str}[{q}];\n",
}
# Get the function from switcher dictionaryresult=model_switcher.get(choice, lambda: "Invalid argument")
print(result)
make_variables(model)
print(f"\tfor (int i = 0; i < {q}; i++) {{\n\t"f"\t{temp_pop_str}[i] = {pop_in_str}[i];}}")
populations=get_DF(q, pop_in_str)
temp_populations=get_DF(q, temp_pop_str)
cm_eq=get_DF(q, cm_eq_pop_str)
F_cm=get_DF(q, F_str)
m=Mraw*temp_populationsprint("\n\t//raw moments from density-probability functions")
# print("\t//[m00, m10, m01, m20, m02, m11, m21, m12, m22]")print_as_vector(m, print_symbol=pop_in_str)
print("\n\t//central moments from raw moments")
cm=Nraw*populationsprint_as_vector(cm, print_symbol=temp_pop_str)
print("\n\t//collision in central moments space")
# print("//calculate equilibrium distributions in cm space")# print("real_t {cm_eq_pop_str}[{q}];\n")# print_as_vector(hardcoded_cm_eq, cm_eq_pop_str) # save time, verbosity# print("//calculate forces in cm space")# print("real_t {F_cm_str}[{q}];")# print_as_vector(hardcoded_F_cm, F_cm_str) # save time, verbosityprint("\t//collide")
defmake_collision(choice):
model_switcher= {
# Relax 2nd moments for hydro, SOI'hydro': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq+ (eye(q) -S_Relax/2) *hardcoded_F_cm,
# Relax 1st moments for ADE, SOI'ade_with_f': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq+ (eye(q) -S_Relax/2) *hardcoded_F_cm,
# Relax 1st moments for ADE, SOI without force'ade': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq,
}
# Get the function from switcher dictionarycm_after_collision=model_switcher.get(choice, lambda: "Invalid argument")
print_as_vector(cm_after_collision, print_symbol=pop_in_str)
make_collision(model)
print("\n\t//back to raw moments")
m=Nraw.inv() *populationsprint_as_vector(m, print_symbol=temp_pop_str)
print("\n\t//back to density-probability functions")
populations=Mraw.inv() *temp_populationsprint_as_vector(populations, print_symbol=pop_in_str)
print("\n}\n")
Output
CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t x_in[9], real_t omega_nu, vector_t u, vector_t F) {
//=== THIS IS AUTOMATICALLY GENERATED CODE ===
real_t uxuy = u.x*u.y;
real_t ux2 = u.x*u.x;
real_t uy2 = u.y*u.y;
real_t m00 = x_in[0] + x_in[1] + x_in[2] + x_in[3] + x_in[4] + x_in[5] + x_in[6] + x_in[7] + x_in[8];
real_t temp[9];for (int i = 0; i < 9; i++) {
temp[i] = x_in[i];}
//raw moments from density-probability functions
x_in[0] = temp[0] + temp[1] + temp[2] + temp[3] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8];
x_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8];
x_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[5] = temp[5] - temp[6] + temp[7] - temp[8];
x_in[6] = temp[5] + temp[6] - temp[7] - temp[8];
x_in[7] = temp[5] - temp[6] - temp[7] + temp[8];
x_in[8] = temp[5] + temp[6] + temp[7] + temp[8];
//central moments from raw moments
temp[0] = x_in[0];
temp[1] = -u.x*x_in[0] + x_in[1];
temp[2] = -u.y*x_in[0] + x_in[2];
temp[3] = ux2*x_in[0] - 2.*u.x*x_in[1] + x_in[3];
temp[4] = uy2*x_in[0] - 2.*u.y*x_in[2] + x_in[4];
temp[5] = uxuy*x_in[0] - u.x*x_in[2] - u.y*x_in[1] + x_in[5];
temp[6] = -ux2*u.y*x_in[0] + ux2*x_in[2] + 2.*uxuy*x_in[1] - 2.*u.x*x_in[5] - u.y*x_in[3] + x_in[6];
temp[7] = -u.x*uy2*x_in[0] + 2.*uxuy*x_in[2] - u.x*x_in[4] + uy2*x_in[1] - 2.*u.y*x_in[5] + x_in[7];
temp[8] = ux2*uy2*x_in[0] - 2.*ux2*u.y*x_in[2] + ux2*x_in[4] - 2.*u.x*uy2*x_in[1] + 4.*uxuy*x_in[5] - 2.*u.x*x_in[7] + uy2*x_in[3] - 2.*u.y*x_in[6] + x_in[8];
//collision in central moments space
//collide
x_in[0] = m00;
x_in[1] = 1/2.*F.x;
x_in[2] = 1/2.*F.y;
x_in[3] = 1/3.*m00*omega_bulk - 1/2.*omega_bulk*temp[3] - 1/2.*omega_bulk*temp[4] - 1/2.*omega_nu*temp[3] + 1/2.*omega_nu*temp[4] + temp[3];
x_in[4] = 1/3.*m00*omega_bulk - 1/2.*omega_bulk*temp[3] - 1/2.*omega_bulk*temp[4] + 1/2.*omega_nu*temp[3] - 1/2.*omega_nu*temp[4] + temp[4];
x_in[5] = -temp[5]*(omega_nu - 1.);
x_in[6] = 1/6.*F.y;
x_in[7] = 1/6.*F.x;
x_in[8] = 1/9.*m00;
//back to raw moments
temp[0] = x_in[0];
temp[1] = u.x*x_in[0] + x_in[1];
temp[2] = u.y*x_in[0] + x_in[2];
temp[3] = ux2*x_in[0] + 2.*u.x*x_in[1] + x_in[3];
temp[4] = uy2*x_in[0] + 2.*u.y*x_in[2] + x_in[4];
temp[5] = uxuy*x_in[0] + u.x*x_in[2] + u.y*x_in[1] + x_in[5];
temp[6] = ux2*u.y*x_in[0] + ux2*x_in[2] + 2.*uxuy*x_in[1] + 2.*u.x*x_in[5] + u.y*x_in[3] + x_in[6];
temp[7] = u.x*uy2*x_in[0] + 2.*uxuy*x_in[2] + u.x*x_in[4] + uy2*x_in[1] + 2.*u.y*x_in[5] + x_in[7];
temp[8] = ux2*uy2*x_in[0] + 2.*ux2*u.y*x_in[2] + ux2*x_in[4] + 2.*u.x*uy2*x_in[1] + 4.*uxuy*x_in[5] + 2.*u.x*x_in[7] + uy2*x_in[3] + 2.*u.y*x_in[6] + x_in[8];
//back to density-probability functions
x_in[0] = temp[0] - temp[3] - temp[4] + temp[8];
x_in[1] = 1/2.*temp[1] + 1/2.*temp[3] - 1/2.*temp[7] - 1/2.*temp[8];
x_in[2] = 1/2.*temp[2] + 1/2.*temp[4] - 1/2.*temp[6] - 1/2.*temp[8];
x_in[3] = -1/2.*temp[1] + 1/2.*temp[3] + 1/2.*temp[7] - 1/2.*temp[8];
x_in[4] = -1/2.*temp[2] + 1/2.*temp[4] + 1/2.*temp[6] - 1/2.*temp[8];
x_in[5] = 1/4.*temp[5] + 1/4.*temp[6] + 1/4.*temp[7] + 1/4.*temp[8];
x_in[6] = -1/4.*temp[5] + 1/4.*temp[6] - 1/4.*temp[7] + 1/4.*temp[8];
x_in[7] = 1/4.*temp[5] - 1/4.*temp[6] - 1/4.*temp[7] + 1/4.*temp[8];
x_in[8] = -1/4.*temp[5] - 1/4.*temp[6] + 1/4.*temp[7] + 1/4.*temp[8];
}
About
Various tools for the TCLB project
Resources
Readme
GPL-3.0 license
Activity
Stars
0 stars
Watchers
1 watching
Forks
0 forks
Report repository
Releases
Packages
Contributors
Languages

, 'i');
if (__m === '*' || __re.test(location.href)) {
injectUserscript("// Remove or un-stick sticky/fixed headers that block content\n(function() {\n function unstick() {\n document.querySelectorAll('header, nav, [role=\"banner\"], .header, .navbar, .sticky, .fixed-top, [style*=\"position: fixed\"], [style*=\"position:sticky\"]').forEach(function(el) {\n if (el.style.position === 'fixed' || el.style.position === 'sticky' || \n getComputedStyle(el).position === 'fixed' || getComputedStyle(el).position === 'sticky') {\n el.style.position = 'static';\n el.style.top = 'auto';\n el.style.zIndex = 'auto';\n }\n });\n }\n \n unstick();\n \n var observer = new MutationObserver(unstick);\n observer.observe(document.body, { childList: true, subtree: true, attributes: true, attributeFilter: ['style', 'class'] });\n})();", "Kill Sticky Headers");
}
} catch(__e) { console.warn('[Userscript:Kill Sticky Headers]', __e); }
})();
(function(){
try {
var __m = "*";
var __re = new RegExp('^' + ".*" + '

Skip to content
BranchesTags
Repository files navigation
TCLB tools
Various tools and extended examples for the TCLB project.
Structure
All tools are divided by language: R, Python, etc. In each directory, the subdirectories should be specific packages related to different aspects of interaction with TCLB
Examples and Papers include scripts used to obtain data used in "real world" cases. Be carefoul - they might be computationally intensive. See Readmi for details
Installation
R
Supervisor: Łukasz Łaniewski-Wołłk
devtools::install_github("CFD-GO/TCLB_tools/R/TCLBtools")
Python
Supervisor: Michał Dzikowski
Content of Python/ directory is meant as a packages repository. Add XXXX/Python to your PYTHONPATH than import packages by subdirectory name.
Python: symbolic_tools
Supervisor: Grzegorz Gruszczyński
It is recomended to use PyCharm to run the scripts.
It facilitates recognition of paths and folders.
Create a PyCharm project by opening the symbolic_tools/ directory in the editor.
Run the examples/ . Some are listed below.
Example: Forcing terms
fromSymbolicCollisions.core.printersimportprint_as_vectorfromSymbolicCollisions.core.ContinousCMTransformsimport \
ContinousCMTransforms, get_mom_vector_from_continuous_deffromSymbolicCollisions.core.cm_symbolsimport \
F3D, dzeta3D, u3DfromSymbolicCollisions.core.cm_symbolsimportrho, moments_dictccmt=ContinousCMTransforms(dzeta3D, u3D, F3D, rho)
lattice='D2Q9'print("\n--- FORCES ---")
print('\n\n// === continuous central moments === \n ')
print('\n//Force -> Force_cm - from continous definition: \n''k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) \n''where fun = forceM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o ')
F_cm=get_mom_vector_from_continuous_def(ccmt.get_force_He_MB,
continuous_transformation=ccmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(F_cm, 'F_cm')
Output:
//Force -> Force_cm - from continous definition: k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) where fun = forceM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o F_cm[0] = 0;
F_cm[1] = Fhydro.x*m00/rho;
F_cm[2] = Fhydro.y*m00/rho;
F_cm[3] = 0;
F_cm[4] = 0;
F_cm[5] = 0;
F_cm[6] = 1/3.*Fhydro.y*m00/rho;
F_cm[7] = 1/3.*Fhydro.x*m00/rho;
F_cm[8] = 0;
Example: Equilibrium Distribution
fromSymbolicCollisions.core.printersimportprint_as_vectorfromsympy.matricesimportMatrixfromsympyimportSymbolfromSymbolicCollisions.core.ContinousCMTransformsimportContinousCMTransforms, get_mom_vector_from_continuous_deffromSymbolicCollisions.core.cm_symbolsimport \
F3D, dzeta3D, u3DfromSymbolicCollisions.core.cm_symbolsimporte_D2Q9, u2D, F2D, rho, moments_dictfromSymbolicCollisions.core.DiscreteCMTransformsimport \
DiscreteCMTransforms, get_mom_vector_from_discrete_deflattice='D2Q9'ccmt=ContinousCMTransforms(dzeta3D, u3D, F3D, rho)
dcmt=DiscreteCMTransforms(e_D2Q9, u2D, F2D, rho)
print("\n--- EQUILIBRIA ---")
print('// === discrete cm ===\n ')
print('\n//population_eq -> cm_eq - by definition: k_mn = sum( (e_ix-ux)^m (e_iy-uy)^n * population_eq_i)')
print("moments: first order (linear) velocity expansion.")
pop_eq=get_mom_vector_from_discrete_def(lambdai: Symbol('m00') *dcmt.get_gamma_first_order(i),
discrete_transform=dcmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(pop_eq, 'cm_eq_first_order')
print("moments: second order (quadratic) velocity expansion.")
pop_eq=get_mom_vector_from_discrete_def(lambdai: Symbol('m00') *dcmt.get_gamma(i),
discrete_transform=dcmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(pop_eq, 'cm_eq_second_order')
print('\n\n// === continous cm === \n ')
# to calculate particular momentrow=moments_dict['D2Q9'][0]
moment=ccmt.get_cm(row, ccmt.get_Maxwellian_DF)
print_as_vector(Matrix([moment]), 'particular_moment')
print('\n//population_eq -> cm_eq - from continous definition: \n''k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) \n''where fun = fM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o ')
cm_eq=get_mom_vector_from_continuous_def(ccmt.get_Maxwellian_DF,
continuous_transformation=ccmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(cm_eq, 'cm_eq')
Output
--- EQUILIBRIA ---
// === discrete cm ===
//population_eq -> cm_eq - by definition: k_mn = sum( (e_ix-ux)^m (e_iy-uy)^n * population_eq_i)
moments: first order (linear) velocity expansion.
cm_eq_first_order[0] = m00;
cm_eq_first_order[1] = 0;
cm_eq_first_order[2] = 0;
cm_eq_first_order[3] = m00*(-ux2 + 1/3.);
cm_eq_first_order[4] = m00*(-uy2 + 1/3.);
cm_eq_first_order[5] = -m00*uxuy;
cm_eq_first_order[6] = 2.*m00*ux2*u.y;
cm_eq_first_order[7] = 2.*m00*u.x*uy2;
cm_eq_first_order[8] = m00*(-3.*ux2*uy2 - 1/3.*ux2 - 1/3.*uy2 + 1/9.);
moments: second order (quadratic) velocity expansion.
cm_eq_second_order[0] = m00;
cm_eq_second_order[1] = 0;
cm_eq_second_order[2] = 0;
cm_eq_second_order[3] = 1/3.*m00;
cm_eq_second_order[4] = 1/3.*m00;
cm_eq_second_order[5] = 0;
cm_eq_second_order[6] = -m00*ux2*u.y;
cm_eq_second_order[7] = -m00*u.x*uy2;
cm_eq_second_order[8] = m00*(3.*ux2*uy2 + 1/9.);
// === continous cm === particular_moment[0] = m00;
//population_eq -> cm_eq - from continous definition: k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) where fun = fM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o cm_eq[0] = m00;
cm_eq[1] = 0;
cm_eq[2] = 0;
cm_eq[3] = 1/3.*m00;
cm_eq[4] = 1/3.*m00;
cm_eq[5] = 0;
cm_eq[6] = 0;
cm_eq[7] = 0;
cm_eq[8] = 1/9.*m00;
Example: Collision Kernel
fromsympy.matricesimporteyefromsympy.printingimportprint_ccodefromSymbolicCollisions.core.cm_symbolsimportomega_ade, omega_b, omega_v, m00fromSymbolicCollisions.core.cm_symbolsimportdynamic_importfromSymbolicCollisions.core.DiscreteCMTransformsimportget_DF, get_m00fromSymbolicCollisions.core.printersimportprint_u2, print_as_vectorfromSymbolicCollisions.core.MatrixGeneratorimportget_raw_moments_matrix, get_shift_matrix# inspired by:# "Consistent Forcing Scheme in the cascaded LBM" L. Fei et al. 2017# eqs 8-12 : (eye(q)-S)*cm + S*cm_eq + (eye(q)-S/2.)*force_in_cm_space# SETUPd=2q=9model='hydro'# choose from '['hydro', 'ade', 'ade_with_f']# DYNAMIC IMPORTSex=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ex_D{d}Q{q}")
ey=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ey_D{d}Q{q}")
ifd==3:
ez=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ez_D{d}Q{q}")
else:
ez=Nonedefget_s_relax_switcher(choice):
s_relax_switcher= {
'hydro': ("SymbolicCollisions.core.cm_symbols", f"S_relax_hydro_D{d}Q{q}"),
'ade_with_f': ("SymbolicCollisions.core.cm_symbols", f"S_relax_ADE_D{d}Q{q}"),
'ade': ("SymbolicCollisions.core.cm_symbols", f"S_relax_ADE_D{d}Q{q}"),
}
which_model=s_relax_switcher.get(choice, lambda: "Invalid argument")
returndynamic_import(*which_model)
S_Relax=get_s_relax_switcher(model)
hardcoded_cm_eq=dynamic_import("SymbolicCollisions.core.hardcoded_results", f"hardcoded_cm_eq_compressible_D{d}Q{q}")
hardcoded_F_cm=dynamic_import("SymbolicCollisions.core.hardcoded_results", f"hardcoded_F_cm_pf_D{d}Q{q}")
fromSymbolicCollisions.core.cm_symbolsimportForce_strasF_str# ARRANGE STUFFMraw=get_raw_moments_matrix(ex, ey, ez)
Nraw=get_shift_matrix(Mraw.inv(), ex, ey, ez)
# from sympy import pprint# pprint(Mraw) # see what you have done# pprint(Nraw)pop_in_str='x_in'# symbol defining populationstemp_pop_str='temp'# symbol defining populationscm_eq_pop_str='cm_eq'# symbol defining populations# GENERATE CODEdefmake_header(choice):
model_switcher= {
'hydro': f"CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t {pop_in_str}[{q}], real_t {omega_v}, vector_t u, vector_t {F_str}) \n{{",
'ade_with_f': f"CudaDeviceFunction void relax_and_collide_ADE_with_F(real_t {pop_in_str}[{q}], real_t {omega_ade}, vector_t u, vector_t {F_str}) \n{{",
'ade': f"CudaDeviceFunction void relax_and_collide_ADE(real_t {pop_in_str}[{q}], real_t {omega_ade}, vector_t u) \n{{",
}
result=model_switcher.get(choice, lambda: "Invalid argument")
print(result)
make_header(model)
print("\t//=== THIS IS AUTOMATICALLY GENERATED CODE ===")
# print(f"real_t {sv} = omega;")# print("real_t bulk_visc = 1./6. ;")# print("real_t {sb} = 1./(3*bulk_visc + 0.5);")# print(f"real_t {sb} = omega_bulk;\n") # s_b = 1./(3*bulk_visc + 0.5)print_u2(d)
print_ccode(get_m00(q, pop_in_str), assign_to=f'\treal_t {m00}')
defmake_variables(choice):
model_switcher= {
'hydro': f"\n\treal_t {temp_pop_str}[{q}];\n",
'ade_with_f': f"\n\treal_t {temp_pop_str}[{q}];\n",
'ade': f"\n\treal_t {temp_pop_str}[{q}];\n",
}
# Get the function from switcher dictionaryresult=model_switcher.get(choice, lambda: "Invalid argument")
print(result)
make_variables(model)
print(f"\tfor (int i = 0; i < {q}; i++) {{\n\t"f"\t{temp_pop_str}[i] = {pop_in_str}[i];}}")
populations=get_DF(q, pop_in_str)
temp_populations=get_DF(q, temp_pop_str)
cm_eq=get_DF(q, cm_eq_pop_str)
F_cm=get_DF(q, F_str)
m=Mraw*temp_populationsprint("\n\t//raw moments from density-probability functions")
# print("\t//[m00, m10, m01, m20, m02, m11, m21, m12, m22]")print_as_vector(m, print_symbol=pop_in_str)
print("\n\t//central moments from raw moments")
cm=Nraw*populationsprint_as_vector(cm, print_symbol=temp_pop_str)
print("\n\t//collision in central moments space")
# print("//calculate equilibrium distributions in cm space")# print("real_t {cm_eq_pop_str}[{q}];\n")# print_as_vector(hardcoded_cm_eq, cm_eq_pop_str) # save time, verbosity# print("//calculate forces in cm space")# print("real_t {F_cm_str}[{q}];")# print_as_vector(hardcoded_F_cm, F_cm_str) # save time, verbosityprint("\t//collide")
defmake_collision(choice):
model_switcher= {
# Relax 2nd moments for hydro, SOI'hydro': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq+ (eye(q) -S_Relax/2) *hardcoded_F_cm,
# Relax 1st moments for ADE, SOI'ade_with_f': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq+ (eye(q) -S_Relax/2) *hardcoded_F_cm,
# Relax 1st moments for ADE, SOI without force'ade': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq,
}
# Get the function from switcher dictionarycm_after_collision=model_switcher.get(choice, lambda: "Invalid argument")
print_as_vector(cm_after_collision, print_symbol=pop_in_str)
make_collision(model)
print("\n\t//back to raw moments")
m=Nraw.inv() *populationsprint_as_vector(m, print_symbol=temp_pop_str)
print("\n\t//back to density-probability functions")
populations=Mraw.inv() *temp_populationsprint_as_vector(populations, print_symbol=pop_in_str)
print("\n}\n")
Output
CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t x_in[9], real_t omega_nu, vector_t u, vector_t F) {
//=== THIS IS AUTOMATICALLY GENERATED CODE ===
real_t uxuy = u.x*u.y;
real_t ux2 = u.x*u.x;
real_t uy2 = u.y*u.y;
real_t m00 = x_in[0] + x_in[1] + x_in[2] + x_in[3] + x_in[4] + x_in[5] + x_in[6] + x_in[7] + x_in[8];
real_t temp[9];for (int i = 0; i < 9; i++) {
temp[i] = x_in[i];}
//raw moments from density-probability functions
x_in[0] = temp[0] + temp[1] + temp[2] + temp[3] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8];
x_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8];
x_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[5] = temp[5] - temp[6] + temp[7] - temp[8];
x_in[6] = temp[5] + temp[6] - temp[7] - temp[8];
x_in[7] = temp[5] - temp[6] - temp[7] + temp[8];
x_in[8] = temp[5] + temp[6] + temp[7] + temp[8];
//central moments from raw moments
temp[0] = x_in[0];
temp[1] = -u.x*x_in[0] + x_in[1];
temp[2] = -u.y*x_in[0] + x_in[2];
temp[3] = ux2*x_in[0] - 2.*u.x*x_in[1] + x_in[3];
temp[4] = uy2*x_in[0] - 2.*u.y*x_in[2] + x_in[4];
temp[5] = uxuy*x_in[0] - u.x*x_in[2] - u.y*x_in[1] + x_in[5];
temp[6] = -ux2*u.y*x_in[0] + ux2*x_in[2] + 2.*uxuy*x_in[1] - 2.*u.x*x_in[5] - u.y*x_in[3] + x_in[6];
temp[7] = -u.x*uy2*x_in[0] + 2.*uxuy*x_in[2] - u.x*x_in[4] + uy2*x_in[1] - 2.*u.y*x_in[5] + x_in[7];
temp[8] = ux2*uy2*x_in[0] - 2.*ux2*u.y*x_in[2] + ux2*x_in[4] - 2.*u.x*uy2*x_in[1] + 4.*uxuy*x_in[5] - 2.*u.x*x_in[7] + uy2*x_in[3] - 2.*u.y*x_in[6] + x_in[8];
//collision in central moments space
//collide
x_in[0] = m00;
x_in[1] = 1/2.*F.x;
x_in[2] = 1/2.*F.y;
x_in[3] = 1/3.*m00*omega_bulk - 1/2.*omega_bulk*temp[3] - 1/2.*omega_bulk*temp[4] - 1/2.*omega_nu*temp[3] + 1/2.*omega_nu*temp[4] + temp[3];
x_in[4] = 1/3.*m00*omega_bulk - 1/2.*omega_bulk*temp[3] - 1/2.*omega_bulk*temp[4] + 1/2.*omega_nu*temp[3] - 1/2.*omega_nu*temp[4] + temp[4];
x_in[5] = -temp[5]*(omega_nu - 1.);
x_in[6] = 1/6.*F.y;
x_in[7] = 1/6.*F.x;
x_in[8] = 1/9.*m00;
//back to raw moments
temp[0] = x_in[0];
temp[1] = u.x*x_in[0] + x_in[1];
temp[2] = u.y*x_in[0] + x_in[2];
temp[3] = ux2*x_in[0] + 2.*u.x*x_in[1] + x_in[3];
temp[4] = uy2*x_in[0] + 2.*u.y*x_in[2] + x_in[4];
temp[5] = uxuy*x_in[0] + u.x*x_in[2] + u.y*x_in[1] + x_in[5];
temp[6] = ux2*u.y*x_in[0] + ux2*x_in[2] + 2.*uxuy*x_in[1] + 2.*u.x*x_in[5] + u.y*x_in[3] + x_in[6];
temp[7] = u.x*uy2*x_in[0] + 2.*uxuy*x_in[2] + u.x*x_in[4] + uy2*x_in[1] + 2.*u.y*x_in[5] + x_in[7];
temp[8] = ux2*uy2*x_in[0] + 2.*ux2*u.y*x_in[2] + ux2*x_in[4] + 2.*u.x*uy2*x_in[1] + 4.*uxuy*x_in[5] + 2.*u.x*x_in[7] + uy2*x_in[3] + 2.*u.y*x_in[6] + x_in[8];
//back to density-probability functions
x_in[0] = temp[0] - temp[3] - temp[4] + temp[8];
x_in[1] = 1/2.*temp[1] + 1/2.*temp[3] - 1/2.*temp[7] - 1/2.*temp[8];
x_in[2] = 1/2.*temp[2] + 1/2.*temp[4] - 1/2.*temp[6] - 1/2.*temp[8];
x_in[3] = -1/2.*temp[1] + 1/2.*temp[3] + 1/2.*temp[7] - 1/2.*temp[8];
x_in[4] = -1/2.*temp[2] + 1/2.*temp[4] + 1/2.*temp[6] - 1/2.*temp[8];
x_in[5] = 1/4.*temp[5] + 1/4.*temp[6] + 1/4.*temp[7] + 1/4.*temp[8];
x_in[6] = -1/4.*temp[5] + 1/4.*temp[6] - 1/4.*temp[7] + 1/4.*temp[8];
x_in[7] = 1/4.*temp[5] - 1/4.*temp[6] - 1/4.*temp[7] + 1/4.*temp[8];
x_in[8] = -1/4.*temp[5] - 1/4.*temp[6] + 1/4.*temp[7] + 1/4.*temp[8];
}
About
Various tools for the TCLB project
Resources
Readme
GPL-3.0 license
Activity
Stars
0 stars
Watchers
1 watching
Forks
0 forks
Report repository
Releases
Packages
Contributors
Languages

, 'i');
if (__m === '*' || __re.test(location.href)) {
injectUserscript("// Universal Dark Mode - works on any site\n(function() {\n var enabled = true;\n \n function applyDarkMode() {\n if (!enabled) return;\n \n // Create style element if it doesn't exist\n var style = document.getElementById('universal-dark-mode-style');\n if (!style) {\n style = document.createElement('style');\n style.id = 'universal-dark-mode-style';\n document.head.appendChild(style);\n }\n \n // Dark mode CSS - inverts colors but preserves images/video\n style.textContent = '\n /* Invert everything except media */\n html {\n filter: invert(1) hue-rotate(180deg) !important;\n background: #1a1a2e !important;\n }\n \n /* Restore images, videos, iframes, canvas */\n img, video, iframe, canvas, svg, picture, [style*=\"background-image\"] {\n filter: invert(1) hue-rotate(180deg) !important;\n }\n \n /* Preserve specific elements that should not be inverted */\n .no-dark-mode, .no-dark-mode *,\n [data-theme=\"light\"], [data-theme=\"light\"],\n .ace_editor, .ace_editor *,\n .CodeMirror, .CodeMirror *,\n .monaco-editor, .monaco-editor *,\n .markdown-body pre, .markdown-body pre *,\n .highlight, .highlight *,\n pre code, pre code * {\n filter: none !important;\n }\n \n /* Fix common UI elements */\n .modal, .popup, .dropdown-menu, .tooltip, .popover {\n filter: invert(1) hue-rotate(180deg) !important;\n background: #2d2d44 !important;\n border-color: #444 !important;\n }\n \n /* Scrollbars */\n ::-webkit-scrollbar { background: #1a1a2e !important; }\n ::-webkit-scrollbar-thumb { background: #444 !important; }\n ::-webkit-scrollbar-thumb:hover { background: #555 !important; }\n \n /* Selection */\n ::selection { background: #4ecdc4 !important; color: #1a1a2e !important; }\n ::-moz-selection { background: #4ecdc4 !important; color: #1a1a2e !important; }\n ';\n }\n \n function removeDarkMode() {\n var style = document.getElementById('universal-dark-mode-style');\n if (style) style.remove();\n }\n \n // Toggle with Alt+Shift+D\n document.addEventListener('keydown', function(e) {\n if (e.altKey && e.shiftKey && e.key === 'D') {\n e.preventDefault();\n enabled = !enabled;\n if (enabled) {\n applyDarkMode();\n console.log('[Universal Dark Mode] Enabled');\n } else {\n removeDarkMode();\n console.log('[Universal Dark Mode] Disabled');\n }\n }\n });\n \n // Apply on load\n applyDarkMode();\n \n // Re-apply on dynamic content\n var observer = new MutationObserver(function(mutations) {\n if (enabled && !document.getElementById('universal-dark-mode-style')) {\n applyDarkMode();\n }\n });\n observer.observe(document.head, { childList: true });\n \n console.log('[Universal Dark Mode] Loaded - Press Alt+Shift+D to toggle');\n})();", "Universal Dark Mode");
}
} catch(__e) { console.warn('[Userscript:Universal Dark Mode]', __e); }
})();
})();

Skip to content
BranchesTags
Repository files navigation
TCLB tools
Various tools and extended examples for the TCLB project.
Structure
All tools are divided by language: R, Python, etc. In each directory, the subdirectories should be specific packages related to different aspects of interaction with TCLB
Examples and Papers include scripts used to obtain data used in "real world" cases. Be carefoul - they might be computationally intensive. See Readmi for details
Installation
R
Supervisor: Łukasz Łaniewski-Wołłk
devtools::install_github("CFD-GO/TCLB_tools/R/TCLBtools")
Python
Supervisor: Michał Dzikowski
Content of Python/ directory is meant as a packages repository. Add XXXX/Python to your PYTHONPATH than import packages by subdirectory name.
Python: symbolic_tools
Supervisor: Grzegorz Gruszczyński
It is recomended to use PyCharm to run the scripts.
It facilitates recognition of paths and folders.
Create a PyCharm project by opening the symbolic_tools/ directory in the editor.
Run the examples/ . Some are listed below.
Example: Forcing terms
fromSymbolicCollisions.core.printersimportprint_as_vectorfromSymbolicCollisions.core.ContinousCMTransformsimport \
ContinousCMTransforms, get_mom_vector_from_continuous_deffromSymbolicCollisions.core.cm_symbolsimport \
F3D, dzeta3D, u3DfromSymbolicCollisions.core.cm_symbolsimportrho, moments_dictccmt=ContinousCMTransforms(dzeta3D, u3D, F3D, rho)
lattice='D2Q9'print("\n--- FORCES ---")
print('\n\n// === continuous central moments === \n ')
print('\n//Force -> Force_cm - from continous definition: \n''k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) \n''where fun = forceM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o ')
F_cm=get_mom_vector_from_continuous_def(ccmt.get_force_He_MB,
continuous_transformation=ccmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(F_cm, 'F_cm')
Output:
//Force -> Force_cm - from continous definition: k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) where fun = forceM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o F_cm[0] = 0;
F_cm[1] = Fhydro.x*m00/rho;
F_cm[2] = Fhydro.y*m00/rho;
F_cm[3] = 0;
F_cm[4] = 0;
F_cm[5] = 0;
F_cm[6] = 1/3.*Fhydro.y*m00/rho;
F_cm[7] = 1/3.*Fhydro.x*m00/rho;
F_cm[8] = 0;
Example: Equilibrium Distribution
fromSymbolicCollisions.core.printersimportprint_as_vectorfromsympy.matricesimportMatrixfromsympyimportSymbolfromSymbolicCollisions.core.ContinousCMTransformsimportContinousCMTransforms, get_mom_vector_from_continuous_deffromSymbolicCollisions.core.cm_symbolsimport \
F3D, dzeta3D, u3DfromSymbolicCollisions.core.cm_symbolsimporte_D2Q9, u2D, F2D, rho, moments_dictfromSymbolicCollisions.core.DiscreteCMTransformsimport \
DiscreteCMTransforms, get_mom_vector_from_discrete_deflattice='D2Q9'ccmt=ContinousCMTransforms(dzeta3D, u3D, F3D, rho)
dcmt=DiscreteCMTransforms(e_D2Q9, u2D, F2D, rho)
print("\n--- EQUILIBRIA ---")
print('// === discrete cm ===\n ')
print('\n//population_eq -> cm_eq - by definition: k_mn = sum( (e_ix-ux)^m (e_iy-uy)^n * population_eq_i)')
print("moments: first order (linear) velocity expansion.")
pop_eq=get_mom_vector_from_discrete_def(lambdai: Symbol('m00') *dcmt.get_gamma_first_order(i),
discrete_transform=dcmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(pop_eq, 'cm_eq_first_order')
print("moments: second order (quadratic) velocity expansion.")
pop_eq=get_mom_vector_from_discrete_def(lambdai: Symbol('m00') *dcmt.get_gamma(i),
discrete_transform=dcmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(pop_eq, 'cm_eq_second_order')
print('\n\n// === continous cm === \n ')
# to calculate particular momentrow=moments_dict['D2Q9'][0]
moment=ccmt.get_cm(row, ccmt.get_Maxwellian_DF)
print_as_vector(Matrix([moment]), 'particular_moment')
print('\n//population_eq -> cm_eq - from continous definition: \n''k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) \n''where fun = fM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o ')
cm_eq=get_mom_vector_from_continuous_def(ccmt.get_Maxwellian_DF,
continuous_transformation=ccmt.get_cm,
moments_order=moments_dict[lattice])
print_as_vector(cm_eq, 'cm_eq')
Output
--- EQUILIBRIA ---
// === discrete cm ===
//population_eq -> cm_eq - by definition: k_mn = sum( (e_ix-ux)^m (e_iy-uy)^n * population_eq_i)
moments: first order (linear) velocity expansion.
cm_eq_first_order[0] = m00;
cm_eq_first_order[1] = 0;
cm_eq_first_order[2] = 0;
cm_eq_first_order[3] = m00*(-ux2 + 1/3.);
cm_eq_first_order[4] = m00*(-uy2 + 1/3.);
cm_eq_first_order[5] = -m00*uxuy;
cm_eq_first_order[6] = 2.*m00*ux2*u.y;
cm_eq_first_order[7] = 2.*m00*u.x*uy2;
cm_eq_first_order[8] = m00*(-3.*ux2*uy2 - 1/3.*ux2 - 1/3.*uy2 + 1/9.);
moments: second order (quadratic) velocity expansion.
cm_eq_second_order[0] = m00;
cm_eq_second_order[1] = 0;
cm_eq_second_order[2] = 0;
cm_eq_second_order[3] = 1/3.*m00;
cm_eq_second_order[4] = 1/3.*m00;
cm_eq_second_order[5] = 0;
cm_eq_second_order[6] = -m00*ux2*u.y;
cm_eq_second_order[7] = -m00*u.x*uy2;
cm_eq_second_order[8] = m00*(3.*ux2*uy2 + 1/9.);
// === continous cm === particular_moment[0] = m00;
//population_eq -> cm_eq - from continous definition: k_mn = integrate(fun, (x, -oo, oo), (y, -oo, oo)) where fun = fM(rho,u,x,y) *(x-ux)^m *(y-uy)^n *(z-uz)^o cm_eq[0] = m00;
cm_eq[1] = 0;
cm_eq[2] = 0;
cm_eq[3] = 1/3.*m00;
cm_eq[4] = 1/3.*m00;
cm_eq[5] = 0;
cm_eq[6] = 0;
cm_eq[7] = 0;
cm_eq[8] = 1/9.*m00;
Example: Collision Kernel
fromsympy.matricesimporteyefromsympy.printingimportprint_ccodefromSymbolicCollisions.core.cm_symbolsimportomega_ade, omega_b, omega_v, m00fromSymbolicCollisions.core.cm_symbolsimportdynamic_importfromSymbolicCollisions.core.DiscreteCMTransformsimportget_DF, get_m00fromSymbolicCollisions.core.printersimportprint_u2, print_as_vectorfromSymbolicCollisions.core.MatrixGeneratorimportget_raw_moments_matrix, get_shift_matrix# inspired by:# "Consistent Forcing Scheme in the cascaded LBM" L. Fei et al. 2017# eqs 8-12 : (eye(q)-S)*cm + S*cm_eq + (eye(q)-S/2.)*force_in_cm_space# SETUPd=2q=9model='hydro'# choose from '['hydro', 'ade', 'ade_with_f']# DYNAMIC IMPORTSex=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ex_D{d}Q{q}")
ey=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ey_D{d}Q{q}")
ifd==3:
ez=dynamic_import("SymbolicCollisions.core.cm_symbols", f"ez_D{d}Q{q}")
else:
ez=Nonedefget_s_relax_switcher(choice):
s_relax_switcher= {
'hydro': ("SymbolicCollisions.core.cm_symbols", f"S_relax_hydro_D{d}Q{q}"),
'ade_with_f': ("SymbolicCollisions.core.cm_symbols", f"S_relax_ADE_D{d}Q{q}"),
'ade': ("SymbolicCollisions.core.cm_symbols", f"S_relax_ADE_D{d}Q{q}"),
}
which_model=s_relax_switcher.get(choice, lambda: "Invalid argument")
returndynamic_import(*which_model)
S_Relax=get_s_relax_switcher(model)
hardcoded_cm_eq=dynamic_import("SymbolicCollisions.core.hardcoded_results", f"hardcoded_cm_eq_compressible_D{d}Q{q}")
hardcoded_F_cm=dynamic_import("SymbolicCollisions.core.hardcoded_results", f"hardcoded_F_cm_pf_D{d}Q{q}")
fromSymbolicCollisions.core.cm_symbolsimportForce_strasF_str# ARRANGE STUFFMraw=get_raw_moments_matrix(ex, ey, ez)
Nraw=get_shift_matrix(Mraw.inv(), ex, ey, ez)
# from sympy import pprint# pprint(Mraw) # see what you have done# pprint(Nraw)pop_in_str='x_in'# symbol defining populationstemp_pop_str='temp'# symbol defining populationscm_eq_pop_str='cm_eq'# symbol defining populations# GENERATE CODEdefmake_header(choice):
model_switcher= {
'hydro': f"CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t {pop_in_str}[{q}], real_t {omega_v}, vector_t u, vector_t {F_str}) \n{{",
'ade_with_f': f"CudaDeviceFunction void relax_and_collide_ADE_with_F(real_t {pop_in_str}[{q}], real_t {omega_ade}, vector_t u, vector_t {F_str}) \n{{",
'ade': f"CudaDeviceFunction void relax_and_collide_ADE(real_t {pop_in_str}[{q}], real_t {omega_ade}, vector_t u) \n{{",
}
result=model_switcher.get(choice, lambda: "Invalid argument")
print(result)
make_header(model)
print("\t//=== THIS IS AUTOMATICALLY GENERATED CODE ===")
# print(f"real_t {sv} = omega;")# print("real_t bulk_visc = 1./6. ;")# print("real_t {sb} = 1./(3*bulk_visc + 0.5);")# print(f"real_t {sb} = omega_bulk;\n") # s_b = 1./(3*bulk_visc + 0.5)print_u2(d)
print_ccode(get_m00(q, pop_in_str), assign_to=f'\treal_t {m00}')
defmake_variables(choice):
model_switcher= {
'hydro': f"\n\treal_t {temp_pop_str}[{q}];\n",
'ade_with_f': f"\n\treal_t {temp_pop_str}[{q}];\n",
'ade': f"\n\treal_t {temp_pop_str}[{q}];\n",
}
# Get the function from switcher dictionaryresult=model_switcher.get(choice, lambda: "Invalid argument")
print(result)
make_variables(model)
print(f"\tfor (int i = 0; i < {q}; i++) {{\n\t"f"\t{temp_pop_str}[i] = {pop_in_str}[i];}}")
populations=get_DF(q, pop_in_str)
temp_populations=get_DF(q, temp_pop_str)
cm_eq=get_DF(q, cm_eq_pop_str)
F_cm=get_DF(q, F_str)
m=Mraw*temp_populationsprint("\n\t//raw moments from density-probability functions")
# print("\t//[m00, m10, m01, m20, m02, m11, m21, m12, m22]")print_as_vector(m, print_symbol=pop_in_str)
print("\n\t//central moments from raw moments")
cm=Nraw*populationsprint_as_vector(cm, print_symbol=temp_pop_str)
print("\n\t//collision in central moments space")
# print("//calculate equilibrium distributions in cm space")# print("real_t {cm_eq_pop_str}[{q}];\n")# print_as_vector(hardcoded_cm_eq, cm_eq_pop_str) # save time, verbosity# print("//calculate forces in cm space")# print("real_t {F_cm_str}[{q}];")# print_as_vector(hardcoded_F_cm, F_cm_str) # save time, verbosityprint("\t//collide")
defmake_collision(choice):
model_switcher= {
# Relax 2nd moments for hydro, SOI'hydro': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq+ (eye(q) -S_Relax/2) *hardcoded_F_cm,
# Relax 1st moments for ADE, SOI'ade_with_f': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq+ (eye(q) -S_Relax/2) *hardcoded_F_cm,
# Relax 1st moments for ADE, SOI without force'ade': (eye(q) -S_Relax) *temp_populations+S_Relax*hardcoded_cm_eq,
}
# Get the function from switcher dictionarycm_after_collision=model_switcher.get(choice, lambda: "Invalid argument")
print_as_vector(cm_after_collision, print_symbol=pop_in_str)
make_collision(model)
print("\n\t//back to raw moments")
m=Nraw.inv() *populationsprint_as_vector(m, print_symbol=temp_pop_str)
print("\n\t//back to density-probability functions")
populations=Mraw.inv() *temp_populationsprint_as_vector(populations, print_symbol=pop_in_str)
print("\n}\n")
Output
CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t x_in[9], real_t omega_nu, vector_t u, vector_t F) {
//=== THIS IS AUTOMATICALLY GENERATED CODE ===
real_t uxuy = u.x*u.y;
real_t ux2 = u.x*u.x;
real_t uy2 = u.y*u.y;
real_t m00 = x_in[0] + x_in[1] + x_in[2] + x_in[3] + x_in[4] + x_in[5] + x_in[6] + x_in[7] + x_in[8];
real_t temp[9];for (int i = 0; i < 9; i++) {
temp[i] = x_in[i];}
//raw moments from density-probability functions
x_in[0] = temp[0] + temp[1] + temp[2] + temp[3] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8];
x_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8];
x_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
x_in[5] = temp[5] - temp[6] + temp[7] - temp[8];
x_in[6] = temp[5] + temp[6] - temp[7] - temp[8];
x_in[7] = temp[5] - temp[6] - temp[7] + temp[8];
x_in[8] = temp[5] + temp[6] + temp[7] + temp[8];
//central moments from raw moments
temp[0] = x_in[0];
temp[1] = -u.x*x_in[0] + x_in[1];
temp[2] = -u.y*x_in[0] + x_in[2];
temp[3] = ux2*x_in[0] - 2.*u.x*x_in[1] + x_in[3];
temp[4] = uy2*x_in[0] - 2.*u.y*x_in[2] + x_in[4];
temp[5] = uxuy*x_in[0] - u.x*x_in[2] - u.y*x_in[1] + x_in[5];
temp[6] = -ux2*u.y*x_in[0] + ux2*x_in[2] + 2.*uxuy*x_in[1] - 2.*u.x*x_in[5] - u.y*x_in[3] + x_in[6];
temp[7] = -u.x*uy2*x_in[0] + 2.*uxuy*x_in[2] - u.x*x_in[4] + uy2*x_in[1] - 2.*u.y*x_in[5] + x_in[7];
temp[8] = ux2*uy2*x_in[0] - 2.*ux2*u.y*x_in[2] + ux2*x_in[4] - 2.*u.x*uy2*x_in[1] + 4.*uxuy*x_in[5] - 2.*u.x*x_in[7] + uy2*x_in[3] - 2.*u.y*x_in[6] + x_in[8];
//collision in central moments space
//collide
x_in[0] = m00;
x_in[1] = 1/2.*F.x;
x_in[2] = 1/2.*F.y;
x_in[3] = 1/3.*m00*omega_bulk - 1/2.*omega_bulk*temp[3] - 1/2.*omega_bulk*temp[4] - 1/2.*omega_nu*temp[3] + 1/2.*omega_nu*temp[4] + temp[3];
x_in[4] = 1/3.*m00*omega_bulk - 1/2.*omega_bulk*temp[3] - 1/2.*omega_bulk*temp[4] + 1/2.*omega_nu*temp[3] - 1/2.*omega_nu*temp[4] + temp[4];
x_in[5] = -temp[5]*(omega_nu - 1.);
x_in[6] = 1/6.*F.y;
x_in[7] = 1/6.*F.x;
x_in[8] = 1/9.*m00;
//back to raw moments
temp[0] = x_in[0];
temp[1] = u.x*x_in[0] + x_in[1];
temp[2] = u.y*x_in[0] + x_in[2];
temp[3] = ux2*x_in[0] + 2.*u.x*x_in[1] + x_in[3];
temp[4] = uy2*x_in[0] + 2.*u.y*x_in[2] + x_in[4];
temp[5] = uxuy*x_in[0] + u.x*x_in[2] + u.y*x_in[1] + x_in[5];
temp[6] = ux2*u.y*x_in[0] + ux2*x_in[2] + 2.*uxuy*x_in[1] + 2.*u.x*x_in[5] + u.y*x_in[3] + x_in[6];
temp[7] = u.x*uy2*x_in[0] + 2.*uxuy*x_in[2] + u.x*x_in[4] + uy2*x_in[1] + 2.*u.y*x_in[5] + x_in[7];
temp[8] = ux2*uy2*x_in[0] + 2.*ux2*u.y*x_in[2] + ux2*x_in[4] + 2.*u.x*uy2*x_in[1] + 4.*uxuy*x_in[5] + 2.*u.x*x_in[7] + uy2*x_in[3] + 2.*u.y*x_in[6] + x_in[8];
//back to density-probability functions
x_in[0] = temp[0] - temp[3] - temp[4] + temp[8];
x_in[1] = 1/2.*temp[1] + 1/2.*temp[3] - 1/2.*temp[7] - 1/2.*temp[8];
x_in[2] = 1/2.*temp[2] + 1/2.*temp[4] - 1/2.*temp[6] - 1/2.*temp[8];
x_in[3] = -1/2.*temp[1] + 1/2.*temp[3] + 1/2.*temp[7] - 1/2.*temp[8];
x_in[4] = -1/2.*temp[2] + 1/2.*temp[4] + 1/2.*temp[6] - 1/2.*temp[8];
x_in[5] = 1/4.*temp[5] + 1/4.*temp[6] + 1/4.*temp[7] + 1/4.*temp[8];
x_in[6] = -1/4.*temp[5] + 1/4.*temp[6] - 1/4.*temp[7] + 1/4.*temp[8];
x_in[7] = 1/4.*temp[5] - 1/4.*temp[6] - 1/4.*temp[7] + 1/4.*temp[8];
x_in[8] = -1/4.*temp[5] - 1/4.*temp[6] + 1/4.*temp[7] + 1/4.*temp[8];
}
About
Various tools for the TCLB project
Resources
Readme
GPL-3.0 license
Activity
Stars
0 stars
Watchers
1 watching
Forks
0 forks
Report repository
Releases
Packages
Contributors
Languages