New parallel version compatible with the first revision of the EFDE/EPDE papwer in CAMWA. (2201f560) · Commits · LBM / lbman

Makefile

+5 −5

Original line number	Diff line number	Diff line
		defaults: main
		defaults: lbman

		CXXFLAGS=-O2 -I. -std=c++11
		LIBS=-lm -lginac

		OBJS = lbm.o core.o main.o coefs.o taylor.o
		OBJS = lbm.o core.o lbman.o taylor.o spatial_pde.o pde_struct2.o defs.o
		DEPS = *.h

		main: $(OBJS) $(DEPS)
		g++ $(CXXFLAGS) $(OBJS) -o main $(LIBS)
		lbman: $(OBJS) $(DEPS)
		g++ $(CXXFLAGS) $(OBJS) -o lbman $(LIBS)

		lbm.o: lbm.cpp defs.h lbm.h
		$(CXX) -MD -MP $(CXXFLAGS2) -c -o $@ $<
		@@ -16,7 +16,7 @@ lbm.o: lbm.cpp defs.h lbm.h
		$(CXX) -MD -MP $(CXXFLAGS) -c -o $@ $<

		clean:
		@rm -rf main
		@rm -rf lbman
		@rm -rf *.o
		@rm -rf *.d

coefs.cpp

deleted100644 → 0

+0 −95

Original line number	Diff line number	Diff line
		#include "coefs.h"

		void COEFS::fill(ex &e)
		{
		if (basis.size()<=0) printf("error: COEFS: no basis %d\n", basis.size());
		else
		{
		SC sc;
		sc.init(basis.size());
		sc.coef=(ex)0;
		sc.base=(ex)1;
		fill_recursive(e, sc, 0); // musts be expanded!!!!
		}
		}

		bool COEFS::fill_recursive(ex &e, SC &sc, int level)
		{
		if (level>basis.size()) return false;
		if (level==basis.size()) // last entry
		{
		ex E=expand(e);
		if (E!=0) // non-zero coefficient
		{
		sc.coef=E;
		v.push_back(sc);
		return true;
		} else
		return false;
		}
		for (int ii=0;ii<order;ii++)
		{
		// hotfix to enforce ordering
		int i = (ii==order-1) ? 0 : ii+1;
		if (sc.order()+i<order)
		{
		ex c = coeff(e, basis[level], i);
		if (c!=0)
		{
		sc.exp[level] = i;
		ex oldbase=sc.base;
		sc.base *= pow(basis[level],i);
		fill_recursive(c, sc, level+1);
		sc.base = oldbase;
		sc.exp[level] = 0;
		}
		}
		}
		return true;
		}

		ex COEFS::check(ex &e)
		{
		// check whether all coefficients were correctly read
		ex expr = 0;
		for (int i=0;i<v.size();i++)
		{
		expr += v[i].coef*v[i].base;
		}
		ex diff=expand(e-expr);
		if (diff!=0)
		{
		cout << "check: nonzero diff = $"<<endl << diff <<"$"<< endl<< endl;
		cout << "source: $"<<endl << e<<endl <<"$"<< endl<< endl;
		cout << "target: $"<<endl << expr<<endl <<"$"<< endl<< endl;
		}
		return diff;
		}

		bool COEFS::compare(SC &sc1, SC &sc2)
		{
		if (sc1.N != sc2.N) return false;
		for (int i=0;i<sc1.N;i++)
		if (sc1.exp[i]!=sc2.exp[i]) return false;
		return true;
		}

		void COEFS::add(COEFS &CF2)
		{
		for (int j=0;j<CF2.v.size();j++)
		{
		// find match for
		int match=-1;
		for (int i=0;i<v.size();i++)
		if (compare(v[i],CF2.v[j])) { match = j; break; }
		if (match<0)
		v.push_back(CF2.v[j]);
		else
		v[match].coef += CF2.v[j].coef;
		}
		}

		void COEFS::copy_basis(COEFS &CF2)
		{
		for (int i=0;i<CF2.basis.size();i++) basis.push_back(CF2.basis[i]);
		}

coefs.h

deleted100644 → 0

+0 −40

Original line number	Diff line number	Diff line
		#pragma once

		#include "defs.h"

		const int MAXSC = 1000; // limit

		struct SC
		{
		int N;
		int exp[MAXSC]; // exponents
		ex coef;
		ex base;
		int order() { int sum=0; for (int i=0;i<N;i++) sum+=exp[i]; return sum; }
		void init(int iN)
		{
		N=iN;
		if (N>=MAXSC) printf("error: MAXSC %d exceeded %d\n",MAXSC, N);
		for (int i=0;i<N;i++) exp[i]=0;
		}
		};


		struct COEFS
		{
		std::vector<SC> v;
		std::vector<ex> basis;
		std::vector<ex> basis_tex;
		std::vector<int> pow_dl;
		std::vector<int> pow_dt;
		int order=6; // default
		void fill(ex &e);
		bool fill_recursive(ex &e, SC &sc, int level);
		ex check(ex &e);
		// adds another COEFS to this
		void add(COEFS &CF2);
		// copy basis
		void copy_basis(COEFS &CF2);
		// compare exp[] for two SCs
		bool compare(SC&sc1, SC&sc2);
		};

core.cpp

+1949 −951

File changed.

Preview size limit exceeded, changes collapsed.

core.h

+101 −44

Original line number	Diff line number	Diff line
		@@ -30,9 +30,11 @@

		#include "defs.h"
		#include "lbm.h"
		#include "coefs.h"
		#include "taylor.h"
		#include "cs.h"
		#include "pde_struct2.h"

		enum { NSE, ADE };

		struct CORE : DEF
		{
		@@ -40,82 +42,137 @@ struct CORE : DEF

		int ORDER=2; // Taylor order
		int collision_model=0;
		int feqorder=-1; // negative = symbolic
		int verbosity=0;
		bool constvel=false;
		bool bigprint;
		bool constmatrixA=true;
		bool numericmatrixA=false;
		bool EFDEprint;
		bool use_tay_vels=true;
		bool use_tay_meqs=true;
		bool use_tay_ms=true;
		bool use_symbolic_B=false;
		bool factor_coefs=true;

		bool use_symbolic_matrices=!true;
		bool use_NSE_strategies=!true; // do not use, the "zeros generator" is much faster!!!
		bool use_NSE_omegas = false; // set omega_0 = omega_1 = ... = omega_d+1 = 1 for NSE to speed up everything
		int pde_strategy = NSE;

		// internal
		bool matrix_A_numeric=false;
		ex matrix_A_coef;



		vct f_mu, f_mueq;
		ex f_rho, f_vel[3];
		lst sub_vels;

		char target_dir[1024];
		char model_desc[1024];

		matrix I, E[MQ], A, B, C, M, Mi, BMi, e[MQ];
		matrix As, dBMi;
		matrix I, E[MQ], A, C, M, Mi, e[MQ];
		symbol *symB;
		matrix *As;
		vector<CS> vecCS;
		int CC[MDIM][MQ];
		matrix *CP;
		ex *trP;
		bool allocated=false;
		bool allocated_tay=false;

		ex *coefs_efde, **coefs_epde;
		ex sym_vels, sym_fs, sym_ms, sym_feqs, **sym_meqs;
		// central values:
		ex sym_m, sym_meq, sym_f, sym_feq;

		ex sym_dl=symbol("\\dl");
		ex sym_dt=symbol("\\dt");

		TAYLOR tay_vels, tay_ms, **tay_meqs;
		symbol **sym_vels;
		TAYLOR **tay_vels;
		symbol sym_vel[MDIM];

		matrix const_feq, const_sub_feq, const_meq;
		matrix feq, meq;
		bool Feq_allocated=false;
		ex sym_cs;
		symbol x=symbol("x"),y=symbol("y"),z=symbol("z"),t=symbol("t");
		symbol dl=symbol("dl","\\dl"), dt=symbol("dt","\\dt");

		ex sym_cs, sym_vel[MDIM];
		// optimization coefs
		symbol symA[MDIM][MQ];
		lst optimal_pars[MQ];
		bool optimal_pars_loaded=false;
		void load_optimal_pars();
		void find_optimal_pars(int pde, int lfrom, int lto);
		void set_optimal_pars(int pde, lst &list);

		void taylor_list_vels(lst &list);
		void taylor_list_fs(lst &list);
		void taylor_list_ms(lst &list);
		void taylor_list_meqs(lst &list);
		void taylor_list_vels_functions(lst &list);
		void taylor_list_vels_moments(lst &list);

		void prepare_EFDE_load_table(symtab &table);

		void desc2str(char*str, int k, int j1, int j2, int j3);
		void desc2str2(char*str, int k, int j1, int j2, int j3);
		void number2str(charstr, const chars1, const char*s2, int number);
		void coord2str(char*str, int x, int y, int z);
		void coord2str(charstr, charstr_tex, int t, int x, int y, int z);
		void desc2str(char s_space_coeff, char s_time_coeff, char* s_space_formula, char* s_time_formula, int l, int i1, int i2, int i3);

		std::string ReplaceString(std::string subject, const std::string& search, const std::string& replace);
		inline void simplify_matrix(matrix &A) { simplify_matrix(A,Q,Q); }
		inline void simplify_matrix(matrix &A, int Nrows, int Ncols) { for (int i=0;i<Nrows;i++) for (int j=0;j<Ncols;j++) A(i,j) = expand(A(i,j)); }
		inline void empty_matrix(matrix &A) { empty_matrix(A,Q,Q); }
		inline void empty_matrix(matrix &A, int Nrows, int Ncols) { for (int i=0;i<Nrows;i++) for (int j=0;j<Ncols;j++) A(i,j)=0;; }

		inline int indO(int a, int b=0, int c=0) { return a + b(ORDER+1) + c(ORDER+1)*(ORDER+1); }
		void Meq_symbolic(matrix &M, int l, int x, int y, int z);
		void Feq_symbolic(matrix &F, int l, int x, int y, int z);
		void Feq_symbolic_Taylor(matrix &F, int l, int x, int y, int z);
		void fill_sym_vels_to_basis(COEFS &CF);

		void define_matrix_A_worker(matrix &_A, ex &u, ex &v, ex &w);
		void define_matrix_A_worker(matrix &_A, ex u, ex v, ex w);
		void define_matrix_A();
		void define_matrix_B();
		void define_matrix_M();
		void define_equilibrium();

		void allocate();
		void allocate_Taylor();

		void log(const char* fmt, ... );

		////////////////////////////////////////////////////////////////////////////////////////////
		//
		// STEP 01: compute LBMAN matrices
		//
		////////////////////////////////////////////////////////////////////////////////////////////

		void load_CP(matrix &P, int l, int i1, int i2, int i3, const char *desc, parser &reader);
		void save_matrix(matrix &P, int l, int i1, int i2, int i3, const char name, const char desc);
		ex load_trP(int l, int i1, int i2, int i3, const char *desc, parser &reader);
		void save_trP(ex trP, int l, int i1, int i2, int i3, const char *desc);


		void compute_coefficients();
		void compute_EFDE_coefs();
		void compute_EPDE_coefs();
		void find_optimal_coefficients(matrix &I_A, ex *pars, int i);
		void compute_LBMAN_B(matrix &B, matrix &X, matrix &BB, matrix dBMi, matrix Xs);

		void compute_LBMAN_matrices_SRT_save();

		bool compute_LBMAN_matrices_symbolic_save_check();
		void compute_LBMAN_matrices_symbolic_save_parallel();
		void compute_LBMAN_matrices_symbolic_load(matrixCP, extrP);
		void compute_LBMAN_matrices_symbolic_load_simplify_and_save();
		bool compute_LBMAN_matrices_symbolic_load_simplify_and_save_check();
		void compute_LBMAN_matrices_load(matrixCP, extrP);

		////////////////////////////////////////////////////////////////////////////////////////////
		//
		// STEP 02: compute EFDE coefficients
		//
		////////////////////////////////////////////////////////////////////////////////////////////
		void compute_EFDE_coefs_save();

		void compute_EFDE_parts_save_parallel();

		bool compute_EFDE_parts_save_parallel_check();

		void compute_EFDE_coefs_save_parallel();
		bool compute_EFDE_coefs_save_parallel_check();

		void compute_EFDE_coefs_symbolic_save();

		void load_EFDE_coefs(matrix &coefs, int pde, int l, int i1, int i2, int i3, const char *desc, parser &reader, bool do_subs);
		void load_EPDE_coefs(matrix &coefs, int i0, int i1, int i2, int i3, const char *desc, parser &reader);

		void export_EFDE_coefs_bigprint(ex &c, ofstream &gout, char* str, bool last);
		void export_EFDE_coefs(int pde);

		////////////////////////////////////////////////////////////////////////////////////////////
		//
		// STEP 03: compute EPDE coefficients
		//
		////////////////////////////////////////////////////////////////////////////////////////////
		void compute_EPDE_coefs_and_save();
		bool compute_EPDE_coefs_and_save_check();

		void export_EPDE_coefs(int pde);
		void export_EPDE_coefs_bigprint(ex &expr, ofstream &fout, char *str);

		CORE(int iDIM, int iQ) : DEF(iDIM, iQ) { sprintf(target_dir,"default"); };
		CORE(int iDIM, int iQ); // : DEF(iDIM, iQ) { sprintf(target_dir,"default"); };
		~CORE();
		};