initial commit

	modified:   README.md
	new file:   core.h
	new file:   core.hpp
	new file:   d1q3.h
	new file:   d2q5.h
	new file:   d2q9.h
	new file:   d3q27.h
	new file:   d3q7.h
	new file:   defs.h
	new file:   defs.hpp
	new file:   export_latex
	new file:   export_latex.cpp
	new file:   export_latex_compare
	new file:   export_latex_compare.cpp
	new file:   lbm.h
	new file:   lbmat
	new file:   lbmat.cpp
	new file:   script_backsubst_searcher.py
	new file:   script_backsubst_searcher_all.sh
	new file:   simplify
	new file:   simplify.cpp
	new file:   tpde.h
	new file:   tpde.hpp

Makefile

+defaults: lbmat simplify export_latex export_latex_compare
+
+CXXFLAGS=-O2 -I. -std=c++11
+LIBS=-lm -lginac 
+
+DEPS = *.h *.hpp
+
+lbmat: $(DEPS) lbmat.cpp
+	$(CXX) $(CXXFLAGS) lbmat.cpp -o lbmat  $(LIBS) 
+
+simplify: $(DEPS) simplify.cpp
+	$(CXX) $(CXXFLAGS) simplify.cpp -o simplify  $(LIBS)
+
+export_latex: $(DEPS) export_latex.cpp
+	$(CXX) $(CXXFLAGS) export_latex.cpp -o export_latex  $(LIBS)
+
+export_latex_compare: $(DEPS) export_latex_compare.cpp
+	$(CXX) $(CXXFLAGS) export_latex_compare.cpp -o export_latex_compare  $(LIBS)
+
+clean:
+	@rm -rf export_latex
+	@rm -rf export_latex_compare
+	@rm -rf lbmat
+	@rm -rf simplify
+	@rm -rf *.o
+	@rm -rf *.d
+

README.md

-# lbmat
+```
+ .----------------.  .----------------.  .----------------.  .----------------.  .----------------. 
+| .--------------. || .--------------. || .--------------. || .--------------. || .--------------. |
+| |   _____      | || |   ______     | || | ____    ____ | || |      __      | || |  _________   | |
+| |  |_   _|     | || |  |_   _ \    | || ||_   \  /   _|| || |     /  \     | || | |  _   _  |  | |
+| |    | |       | || |    | |_) |   | || |  |   \/   |  | || |    / /\ \    | || | |_/ | | \_|  | |
+| |    | |   _   | || |    |  __'.   | || |  | |\  /| |  | || |   / ____ \   | || |     | |      | |
+| |   _| |__/ |  | || |   _| |__) |  | || | _| |_\/_| |_ | || | _/ /    \ \_ | || |    _| |_     | |
+| |  |________|  | || |  |_______/   | || ||_____||_____|| || ||____|  |____|| || |   |_____|    | |
+| |              | || |              | || |              | || |              | || |              | |
+| '--------------' || '--------------' || '--------------' || '--------------' || '--------------' |
+ '----------------'  '----------------'  '----------------'  '----------------'  '----------------' 
+```
+
+# Lattice Boltzmann Method Analysis Tool
+
+ https://mmg-gitlab.fjfi.cvut.cz/gitlab/lbm/lbmat
+
+ Developped by:
+* Radek Fučík (fucik@fjfi.cvut.cz) FNSPE CTU in Prague, Czech Republic
+* Robert Straka (straka@metal.agh.edu.pl) AGH Krakow, Poland
+* Jakub Klinkovský (klinkovsky@mmg.fjfi.cvut.cz) FNSPE CTU in Prague, Czech Republic
+* Pavel Eichler (eichlpa1@fjfi.cvut.cz) FNSPE CTU in Prague, Czech Republic
+* Tomáš Oberhuber (oberhuber.tomas@fjfi.cvut.cz) FNSPE CTU in Prague, Czech Republic
+
+ Read and cite: <br/>
+* Radek Fučík, Pavel Eichler, Jakub Klinkovský, Robert Straka and Tomáš Oberhuber: <i>Lattice Boltzmann Method Analysis Tool (LBMAT)</i> in review in Elsevier Computers & Mathematics with Applications
+
+* Radek Fučík and Robert Straka: <i>Equivalent finite difference and partial differential equations for the lattice Boltzmann method</i>, Elsevier Computers & Mathematics with Applications, vol. 90, 2021, 96-103 
+ https://doi.org/10.1016/j.camwa.2021.03.014
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+ https://www.gnu.org/licenses/
 

core.h

+// Lattice Boltzmann Method Analysis Tool
+// https://mmg-gitlab.fjfi.cvut.cz/gitlab/lbm/lbmat
+
+#pragma once
+
+#include "defs.h"
+#include "lbm.h"
+#include "tpde.h"
+
+template< typename LBM_TYPE >
+struct CORE : BASE
+{
+	using LBM=LBM_TYPE;
+	LBM lbm;
+
+	const int MAXO=5;	// maximal time order to deal with, increase this if you need more than that
+	int ORDER=4;		// Taylor order
+	int MFCO=4;		// maximal factors common order for TPDE 
+	int verbosity=0;
+	// link to lbm
+	int Q=0, D=0;
+	bool use_symbolic_coefs=false; // use symbolic coefficients for abs_rhs to save memory
+
+	FUN fun;
+
+	int dump_lst_count=0;
+	int subs_lst_level=0;
+
+	// virtualize for non-matrix LBMs!!! (such as CuLBM...)
+	virtual void assemble_pdes(TPDE *pde);
+
+	bool find_lowest_time_derivative(TPDE &expression, int basindex, int order, int& kt, int& kx, int& ky, int&kz);
+
+	void subs_time_derivatives_up_to_order(TPDE &expr, int basindex, TPDE *ders, int max_order);
+	void eliminate_all_spatial_orders_for_given_time_order(int basindex, int time_order, TPDE *ders);
+	void eliminate_all_spatial_orders_for_given_time_order(TPDE &target, int basindex, int time_order, TPDE &source);
+	void replace_time_derivative(TPDE &target, int bj, int time_order, TPDE &derivative);
+
+	void eliminate_higher_raw_moments(vct& bas, TPDE *pde, TPDE*abs_rhs, int offset);
+
+	// Algorithm 1
+	void gauss_eliminate_derivatives_on_diagonal(int bi, vct& bas, TPDE*abs_rhs, int index_offset);
+	void gauss_eliminate_derivatives_below_diagonal(int bi, vct& bas, TPDE*abs_rhs, int index_offset);
+	void gauss_eliminate_derivatives_above_diagonal(int bi, vct& bas, TPDE*abs_rhs, int index_offset);
+	void gauss_eliminate_coefs_on_diagonal(int bi, vct& bas, TPDE*abs_rhs, int index_offset);
+	void gauss_eliminate_coefs_below_diagonal(int bi, vct& bas, TPDE*abs_rhs, int index_offset);
+	void gauss_eliminate_coefs_above_diagonal(int bi, vct& bas, TPDE*abs_rhs, int index_offset);
+
+	// Algorithm 2
+	void compute_NSE();
+	void compute_ADE();
+
+	void compute_all();
+
+	CORE();
+};
+
+#include "core.hpp"

d1q3.h

+#pragma once
+#include "lbm.h"
+
+template< typename MM>
+struct D1Q3 : LBM<MM>
+{
+	static constexpr const char* id = "d1q3";
+	void setC(matrix &C)
+	{
+		C = matrix(1, 3);
+		int maxima_matrix[1][3] = {
+			{0,	1,	-1}
+		};
+		for (int i=0;i<1;i++)
+		for (int j=0;j<3;j++)
+			C(i,j) = (ex)(maxima_matrix[i][j]);
+	}
+	D1Q3() : LBM<MM>(3) {}
+};
+
+struct D1Q3_M1
+{
+	static constexpr const char* id = "M1";
+	static void setM(matrix &M, matrix &C)
+	{
+		const int Q=3,D=1;
+		M = matrix(Q, Q);
+		for (int j=0;j<Q;j++)
+		{
+			M(0,j) = (ex)1;
+			for (int d=0;d<D;d++) M(d+1,j) = C(d,j);
+			M(2,j) = pow(C(0,j), 2);
+		}
+	}
+	static ex get_moments_eq(int q, FUN& fun, ex sym_cs, matrix &C)
+	{
+		if (q==0) return fun.rho;
+		if (q==1) return fun.rho*fun.vel[0];
+		if (q==2) return fun.rho*(fun.vel[0]*fun.vel[0] + sym_cs*sym_cs);
+		return 0;
+	}
+};
+
+template< typename MM>
+struct D1Q3_SRT : D1Q3<MM>
+{
+	static constexpr const char* id_col = "srt";
+	using LBM<MM>::Q;
+
+	virtual void setA(matrix &A, ex& u, ex&v, ex&w)
+	{
+		A = matrix(Q, Q);
+		for (int i=0;i<Q;i++)
+		for (int j=0;j<Q;j++)
+			A(i,j) = (i==j) ? (ex)LBM<MM>::sym_o : (ex)0;
+	}
+
+};
+
+template< typename MM>
+struct D1Q3_MRT1 : D1Q3<MM>
+{
+	static constexpr const char* id_col = "mrt1";
+	virtual void setM(matrix &M, matrix &C)
+	{
+		D1Q3_M1::setM(M,C);
+	}
+	virtual void setA(matrix &A, ex& u, ex&v, ex&w)
+	{
+		matrix S,M;
+		LBM<MM>::setS(S);
+	        LBM<MM>::setM(M);
+	        matrix Mi = inverse(M);
+	        A = Mi.mul(S).mul(M);
+	}
+};
+
+template< typename MM>
+struct D1Q3_CLBM1 : D1Q3_MRT1<MM>
+{
+	static constexpr const char* id_col = "clbm1";
+	virtual void setA(matrix &A, ex& u, ex&v, ex&w)
+	{
+		matrix S,K;
+		LBM<MM>::setS(S);
+		LBM<MM>::setK(K, u, v, w);
+		matrix Ki = inverse(K);
+		A = Ki.mul(S).mul(K);
+	}
+};

d2q5.h

+#pragma once
+#include "lbm.h"
+
+template< typename MM>
+struct D2Q5 : LBM<MM>
+{
+	static constexpr const char* id = "d2q5";
+
+	virtual void setC(matrix &C)
+	{
+		C = matrix(2, 5);
+		int maxima_matrix[2][5] = {
+			{0,	1,	0,	-1,	0},
+			{0,	0,	1,	0,	-1}
+		};
+		for (int i=0;i<2;i++)
+		for (int j=0;j<5;j++)
+			C(i,j) = (ex)(maxima_matrix[i][j]);
+	}
+	D2Q5() : LBM<MM>(5) {}
+};
+
+struct D2Q5_M1
+{
+	static constexpr const char* id = "M1";
+	static void setM(matrix &M, matrix &C)
+	{
+		const int Q=5,D=2;
+		M = matrix(Q, Q);
+		for (int j=0;j<Q;j++)
+		{
+			M(0,j) = (ex)1;
+			for (int d=0;d<D;d++) M(d+1,j) = C(d,j);
+			M(3,j) = pow(C(0,j), 2);
+			M(4,j) = pow(C(1,j), 2);
+		}
+	}
+	static ex get_moments_eq(int q, FUN& fun, ex sym_cs, matrix &C)
+	{
+		// equilibrium Asinari_PRE_2008bR.pdf P7 Eq 16
+
+		if (q==0) return fun.rho;
+		// beware of correct order_eq of moments - Asinari has different 4th, 5th and 6th row!
+		if (q==1) return fun.rho*fun.vel[0];
+		if (q==2) return fun.rho*fun.vel[1];
+ 		if (q==3) return fun.rho*(fun.vel[0]*fun.vel[0] + sym_cs*sym_cs);
+ 		if (q==4) return fun.rho*(fun.vel[1]*fun.vel[1] + sym_cs*sym_cs);
+		return 0;
+	}
+
+};
+
+struct D2Q5_M2
+{
+	static constexpr const char* id = "M2";
+	static void setM(matrix &M, matrix &C)
+	{
+		const int Q=5,D=2;
+		M = matrix(Q, Q);
+		for (int j=0;j<Q;j++)
+		{
+			// zeroth moment
+			M(0,j) = (ex)1;
+			// first moments
+			for (int d=0;d<D;d++) M(d+1,j) = C(d,j);
+			M(3,j) = pow(C(0,j), 2) + pow(C(1,j), 2);
+			M(4,j) = pow(C(0,j), 2) - pow(C(1,j), 2);
+		}
+	}
+
+	static ex get_moments_eq(int q, FUN&fun, ex sym_cs, matrix &C)
+	{
+		const int Q=5;
+		matrix mu1_eq = matrix(Q,1);
+		for (int i=0;i<Q;i++)
+			mu1_eq(i,0) = D2Q5_M1::get_moments_eq(i,fun,sym_cs,C);
+		matrix M1,M2;
+		D2Q5_M1::setM(M1,C);
+		setM(M2,C);
+		matrix M1i = inverse(M1);
+		matrix mu2_eq = M2.mul(M1i).mul(mu1_eq);
+		return mu2_eq(q,0);
+	}
+};
+
+template< typename MM>
+struct D2Q5_SRT : D2Q5<MM>
+{
+	static constexpr const char* id_col = "srt";
+	using LBM<MM>::Q;
+
+	virtual void setA(matrix &A, ex& u, ex&v, ex&w)
+	{
+		A = matrix(Q, Q);
+		for (int i=0;i<Q;i++)
+		for (int j=0;j<Q;j++)
+			A(i,j) = (i==j) ? (ex)LBM<MM>::sym_o : (ex)0;
+	}
+
+};
+
+template< typename MM>
+struct D2Q5_MRT1 : D2Q5<MM>
+{
+	static constexpr const char* id_col = "mrt1";
+
+	virtual void setM(matrix &M, matrix &C)
+	{
+		D2Q5_M1::setM(M,C);
+	}
+
+	virtual void setA(matrix &A, ex& u, ex&v, ex&w)
+	{
+		matrix S,M;
+		LBM<MM>::setS(S);
+	        LBM<MM>::setM(M);
+	        matrix Mi = inverse(M);
+	        A = Mi.mul(S).mul(M);
+	}
+};
+
+template< typename MM>
+struct D2Q5_MRT2 : D2Q5_MRT1<MM>
+{
+	static constexpr const char* id_col = "mrt2";
+
+	virtual void setM(matrix &M, matrix &C)
+	{
+		D2Q5_M2::setM(M,C);
+	}
+};
+
+template< typename MM>
+struct D2Q5_CLBM1 : D2Q5_MRT1<MM>
+{
+	static constexpr const char* id_col = "clbm1";
+
+	virtual void setA(matrix &A, ex& u, ex&v, ex&w)
+	{
+		matrix S,K;
+		LBM<MM>::setS(S);
+		LBM<MM>::setK(K, u, v, w);
+		matrix Ki = inverse(K);
+		A = Ki.mul(S).mul(K);
+	}
+};
+
+template< typename MM>
+struct D2Q5_CLBM2 : D2Q5_CLBM1<MM>
+{
+	static constexpr const char* id_col = "clbm2";
+
+	virtual void setM(matrix &M, matrix &C)
+	{
+		D2Q5_M2::setM(M,C);
+	}
+};

d2q9.h

+#pragma once
+#include "lbm.h"
+
+template< typename MM>
+struct D2Q9 : LBM<MM>
+{
+	static constexpr const char* id = "d2q9";
+	virtual void setC(matrix &C)
+	{
+		C = matrix(2, 9);
+		int maxima_matrix[2][9] = {
+			{0,	1,	0,	-1,	0,	1,	-1,	-1,	1},
+			{0,	0,	1,	0,	-1,	1,	1,	-1,	-1}
+		};
+		for (int i=0;i<2;i++)
+		for (int j=0;j<9;j++)
+			C(i,j) = (ex)(maxima_matrix[i][j]);
+	}
+	D2Q9() : LBM<MM>(9) {}
+};
+
+struct D2Q9_M1
+{
+	static constexpr const char* id = "M1";
+	static void setM(matrix &M, matrix &C)
+	{
+		const int Q=9,D=2;
+		M = matrix(Q, Q);
+		for (int j=0;j<Q;j++)
+		{
+			// zeroth moment
+			 M(0,j) = (ex)1;
+			// first moments
+			for (int d=0;d<D;d++) M(d+1,j) = C(d,j);
+			// the rest of moments
+			M(3,j) = C(0,j)*C(1,j);
+			M(4,j) = pow(C(0,j), 2);
+			M(5,j) = pow(C(1,j), 2);
+			M(6,j) = pow(C(0,j), 2)*C(1,j);
+			M(7,j) = pow(C(1,j), 2)*C(0,j);
+			M(8,j) = pow(C(0,j), 2)*pow(C(1,j), 2);
+		}
+	}
+
+	static ex get_moments_eq(int q, FUN& fun, ex sym_cs, matrix &C)
+	{
+		// equilibrium Asinari_PRE_2008bR.pdf P7 Eq 16
+		if (q==0) return fun.rho;
+		// beware of correct order_eq of moments - Asinari has different 4th, 5th and 6th row!
+		if (q==1) return fun.rho*fun.vel[0];
+		if (q==2) return fun.rho*fun.vel[1];
+ 		if (q==3) return fun.rho*(fun.vel[0]*fun.vel[1]);
+ 		if (q==4) return fun.rho*(fun.vel[0]*fun.vel[0] + sym_cs*sym_cs);
+ 		if (q==5) return fun.rho*(fun.vel[1]*fun.vel[1] + sym_cs*sym_cs);
+ 		if (q==6) return fun.rho*(fun.vel[1]*sym_cs*sym_cs + fun.vel[0]*fun.vel[0]*fun.vel[1]);
+ 		if (q==7) return fun.rho*(fun.vel[0]*sym_cs*sym_cs + fun.vel[0]*fun.vel[1]*fun.vel[1]);
+		if (q==8) return fun.rho*(sym_cs*sym_cs*sym_cs*sym_cs + sym_cs*sym_cs*(fun.vel[0]*fun.vel[0] + fun.vel[1]*fun.vel[1]) + fun.vel[0]*fun.vel[0]*fun.vel[1]*fun.vel[1]);
+		return 0;
+	}
+
+};
+
+struct D2Q9_M2
+{
+	static constexpr const char* id = "M2";
+	static void setM(matrix &M, matrix &C)
+	{
+		const int Q=9,D=2;
+		M = matrix(Q, Q);
+		for (int j=0;j<Q;j++)
+		{
+			// zeroth moment
+			 M(0,j) = (ex)1;
+			// first moments
+			for (int d=0;d<D;d++) M(d+1,j) = C(d,j);
+			// dev to see whether we can get the same result as in maxima: 
+			// central moment basis k00,k10,k01,k10*k01,k20+k02,k20-k02,k21,k12,k22
+			M(3,j) = C(0,j)*C(1,j);
+			M(4,j) = pow(C(0,j), 2) + pow(C(1,j), 2);
+			M(5,j) = pow(C(0,j), 2) - pow(C(1,j), 2);
+			M(6,j) = pow(C(0,j), 2)*C(1,j);
+			M(7,j) = pow(C(1,j), 2)*C(0,j);
+			M(8,j) = pow(C(0,j), 2)*pow(C(1,j), 2);
+		}
+	}
+
+	static ex get_moments_eq(int q, FUN&fun, ex sym_cs, matrix &C)
+	{
+		const int Q=9;
+		// workaround:
+		// mu_eq = 
+		matrix mu1_eq = matrix(Q,1);
+		for (int i=0;i<Q;i++)
+			mu1_eq(i,0) = D2Q9_M1::get_moments_eq(i,fun,sym_cs,C);
+		matrix M1,M2;
+		D2Q9_M1::setM(M1,C);
+		setM(M2,C);
+		matrix M1i = inverse(M1);
+		matrix mu2_eq = M2.mul(M1i).mul(mu1_eq);
+		return mu2_eq(q,0);
+	}
+};
+
+struct D2Q9_M3dev
+{
+	static constexpr const char* id = "M3dev";
+	static void setM(matrix &M, matrix &C)
+	{
+		const int Q=9,D=2;
+		M = matrix(Q, Q);
+		for (int j=0;j<Q;j++)
+		{
+			// zeroth moment
+			 M(0,j) = (ex)1;
+			// first moments
+			for (int d=0;d<D;d++) M(d+1,j) = C(d,j);
+			for (int d=D+1;d<Q;d++) M(d,j) = 0;
+			if (j>D) M(j,j) = (ex)1;
+		}
+	}
+
+	static ex get_moments_eq(int q, FUN&fun, ex sym_cs, matrix &C)
+	{
+		const int Q=9;
+		// workaround:
+		// mu_eq = 
+		matrix mu1_eq = matrix(Q,1);
+		for (int i=0;i<Q;i++)
+			mu1_eq(i,0) = D2Q9_M1::get_moments_eq(i,fun,sym_cs,C);
+		matrix M1,M2;
+		D2Q9_M1::setM(M1,C);
+		setM(M2,C);
+		matrix M1i = inverse(M1);
+		matrix mu2_eq = M2.mul(M1i).mul(mu1_eq);
+		return mu2_eq(q,0);
+	}
+};
+
+template< typename MM>
+struct D2Q9_SRT : D2Q9<MM>
+{
+	static constexpr const char* id_col = "srt";
+	using LBM<MM>::Q;
+
+	virtual void setA(matrix &A, ex& u, ex&v, ex&w)
+	{
+		A = matrix(Q, Q);
+		for (int i=0;i<Q;i++)
+		for (int j=0;j<Q;j++)
+			A(i,j) = (i==j) ? (ex)LBM<MM>::sym_o : (ex)0;
+	}
+
+};
+
+template< typename MM>
+struct D2Q9_MRT1 : D2Q9<MM>
+{
+	static constexpr const char* id_col = "mrt1";
+	virtual void setM(matrix &M, matrix &C)
+	{
+		D2Q9_M1::setM(M,C);
+	}
+	virtual void setA(matrix &A, ex& u, ex&v, ex&w)
+	{
+		matrix S,M;
+		LBM<MM>::setS(S);
+	        LBM<MM>::setM(M);
+	        matrix Mi = inverse(M);
+	        A = Mi.mul(S).mul(M);
+	}
+};
+
+template< typename MM>
+struct D2Q9_MRT2 : D2Q9_MRT1<MM>
+{
+	static constexpr const char* id_col = "mrt2";
+	virtual void setM(matrix &M, matrix &C)
+	{
+		D2Q9_M2::setM(M,C);
+	}
+};
+
+template< typename MM>
+struct D2Q9_CLBM1 : D2Q9_MRT1<MM>
+{
+	static constexpr const char* id_col = "clbm1";
+	virtual void setA(matrix &A, ex& u, ex&v, ex&w)
+	{
+		matrix S,K;
+		LBM<MM>::setS(S);
+		LBM<MM>::setK(K, u, v, w);
+		matrix Ki = inverse(K);
+		A = Ki.mul(S).mul(K);
+	}
+};
+
+template< typename MM>
+struct D2Q9_CLBM2 : D2Q9_CLBM1<MM>
+{
+	static constexpr const char* id_col = "clbm2";
+	virtual void setM(matrix &M, matrix &C)
+	{
+		D2Q9_M2::setM(M,C);
+	}
+};
+
+template< typename MM>
+struct D2Q9_CuLBM1 : D2Q9<MM>
+{
+    static constexpr const char* id_col = "culbm1";
+    static constexpr const bool is_matrix_lbm = false;
+    using LBM<MM>::sym_os;
+    using LBM<MM>::sym_cs;
+    virtual void collision(matrix &rhs, matrix &dfs, ex& u, ex&v, ex&w, ex&rho)
+    {
+	// converted from lbm: lbm_col_cum.h version 2021.08.21
+	// simplified way of using fractions in the CuLBM codes
+	const ex no1 = (ex)1;
+	const ex n2o3 = 2/(ex)3;
+	const ex n1o2 = 1/(ex)2;