Skip to content

Commit 633612c

Browse files
GalacticBobsterGalacticBobster
committed
Added tridiagonal solver
1 parent b35e21a commit 633612c

File tree

2 files changed

+115
-51
lines changed

2 files changed

+115
-51
lines changed

‎src/tridiagonal_solver.cpp‎

Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
// output stream
2+
#include <iostream>
3+
#include <fstream>
4+
#include <cantera/numerics/eigen_dense.h>
5+
#include <cantera/numerics/eigen_sparse.h>
6+
#include <string>
7+
8+
//We use the Thomas algorithm for solving the block tridiagonal matrix system
9+
Eigen::MatrixXd ThomasSolver(Eigen::MatrixXd A, Eigen::MatrixXd B, int nsp, int nSize){
10+
std::cout << "Starting Thomas solver" << std::endl;
11+
//Initialize output vector
12+
Eigen::MatrixXd Xvec = Eigen::VectorXd::Zero(nsp*nSize);
13+
std::vector<Eigen::MatrixXd> C_prime(nSize);
14+
std::vector<Eigen::MatrixXd> D_prime(nSize);
15+
std::cout << "Initialization done!" << std::endl;
16+
17+
for(int i = 0; i < nSize; i++){
18+
19+
20+
//Forward elimination
21+
if(i == 0){
22+
Eigen::MatrixXd B1 = A.block(0, 0, nsp, nsp);
23+
std::cout << "B1 done" << std::endl;
24+
Eigen::MatrixXd C1 = A.block(0, nsp, nsp, nsp);
25+
std::cout << "C1 done" << std::endl;
26+
Eigen::MatrixXd D1 = B.block(0, 0, nsp, 1);
27+
std::cout << "D1 done" << std::endl;
28+
29+
C_prime[0] = B1.inverse()*C1;
30+
D_prime[0] = B1.inverse()*D1;
31+
std::cout << "Forward elimination at point 0 complete" << std::endl;
32+
}
33+
34+
if(i > 0 && i < nSize - 1){
35+
Eigen::MatrixXd Bi = A.block(i*nsp, i*nsp, nsp, nsp);
36+
std::cout << "Bi done" << std::endl;
37+
Eigen::MatrixXd Ai = A.block((i-1)*nsp, i*nsp, nsp, nsp);
38+
std::cout << "Ai done" << std::endl;
39+
Eigen::MatrixXd Ci = A.block(i*nsp, (i+1)*nsp, nsp, nsp);
40+
std::cout << "Ci done" << std::endl;
41+
Eigen::MatrixXd Di = B.block(i*nsp, 0, nsp, 1);
42+
std::cout << "Di done" << std::endl;
43+
44+
C_prime[i] = (Bi - Ai*C_prime[i-1]).inverse()*Ci;
45+
std::cout << "C prime done" << std::endl;
46+
D_prime[i] = (Bi - Ai*C_prime[i-1]).inverse()*(Di - Ai*D_prime[i-1]);
47+
std::cout << "D prime done" << std::endl;
48+
49+
std::cout << "Forward elimination at point " << i << " complete" << std::endl;
50+
51+
}
52+
53+
if(i == nSize - 1){
54+
Eigen::MatrixXd Bi = A.block(i*nsp, i*nsp, nsp, nsp);
55+
std::cout << "Bi done" << std::endl;
56+
Eigen::MatrixXd Ai = A.block((i-1)*nsp, i*nsp, nsp, nsp);
57+
std::cout << "Ai done" << std::endl;
58+
Eigen::MatrixXd Di = B.block(i*nsp, 0, nsp, 1);
59+
std::cout << "Di done" << std::endl;
60+
61+
D_prime[i] = (Bi - Ai*C_prime[i-1]).inverse()*(Di - Ai*D_prime[i-1]);
62+
std::cout << "D prime done" << std::endl;
63+
64+
std::cout << "Forward elimination at point " << i << " complete" << std::endl;
65+
66+
}
67+
68+
69+
}
70+
std::cout << "Forward elimination complete!" << std::endl;
71+
//Back substitution
72+
std::cout << "Starting backward substitution" << std::endl;
73+
74+
for(int i = nSize-1; i >= 0; i--){
75+
std::cout << "D_prime: " << D_prime[i] << std::endl;
76+
if(i == nSize-1){
77+
Xvec.block(i*nsp, 0, nsp, 1) = D_prime[i] ;
78+
//std::cout << D_prime[i] << std::endl;
79+
std::cout << "Xvec complete" << std::endl;
80+
}
81+
82+
else{
83+
84+
Eigen::MatrixXd X_next = Xvec.block((i+1)*nsp, 0, nsp, 1);
85+
Xvec.block(i*nsp, 0, nsp, 1) = D_prime[i] - C_prime[i]*X_next;
86+
std::cout << Xvec.block(i*nsp, 0, nsp, 1) << std::endl;
87+
88+
}
89+
90+
}
91+
92+
93+
94+
return Xvec;
95+
}

‎tests/1DPP.cpp‎

Lines changed: 20 additions & 51 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,7 @@
3131
#include <interpolation.hpp>
3232
#include <CustomRate.hpp>
3333
#include <CustomTransport.hpp>
34+
#include <tridiagonal_solver.cpp>
3435

3536
// NetCDF Output
3637
#if NETCDFOUTPUT
@@ -625,12 +626,10 @@ while(Ttot < Tmax) {
625626

626627
//The transmission coefficient from opacity
627628
Opacity.col(j) = Opacity.col(j) - (dh*handleCustomOpacity(PlanetName, gasThermo, Press, Temp, AtmData(iAlt,j), stellar_input.row(0)));
628-
//OpacityStorage.col(j) = Opacity.col(j);
629629
Opacity.col(j) = Opacity.col(j).array().exp().matrix();
630630
Opacity.col(j) = (Opacity.col(j).array()*Opacity.col(j-1).array()).matrix();
631631
//Stellar spectrum at each altitude
632632
Stellar_activity.col(j) = (Stellar_activity.col(j-1).array()*Opacity.col(j).array()).transpose().matrix();
633-
// std::cout << stellar_activity.col(j) << std::endl;
634633
}
635634
if(j == 0){
636635
clden_species_list = pinput->GetOrAddString("coldensity", "species", "nan");
@@ -645,14 +644,11 @@ while(Ttot < Tmax) {
645644
species_inx = gas2->speciesIndex(x);
646645
double number_density = Press/(Kb*Temp);
647646
Opacity.col(j) = Opacity.col(j) - (absorber_cross_data.col(Absorber)*cldMag_diff(species_inx)/cos(sz_angle*3.14/180));
648-
//std::cout << absorber_cross_data.col(Absorber).transpose()*cldMag_diff(species_inx)/cos(sz_angle*3.14/180) << std::endl;
649-
// std::cout << Opacity.col(j).transpose() << std::endl;
650647
Absorber++;
651648
}
652649
absorber_species_list = m.suffix().str();
653650
}
654651

655-
//std::cout << "Absorption done! " << std::endl;
656652
//Rayleigh scattering
657653
int Scatter = 0;
658654
std::string scat_species_list = pinput->GetString("scatcross", "scatterers");
@@ -665,29 +661,23 @@ while(Ttot < Tmax) {
665661
}
666662
scat_species_list = m.suffix().str();
667663
}
668-
// std::cout << "Scattering done at TOA!" << std::endl;
669664
//Updating the stellar activity at upper boundary
670665
//The transmission coefficient from opacity
671-
// std::cout << Opacity.col(j) << std::endl;
672-
// OpacityStorage.col(j) = Opacity.col(j);
673666
Opacity.col(j) = Opacity.col(j).array().exp().matrix();
674667
//Stellar spectrum at each altitude
675668
Stellar_activity.col(j) = (stellar_input.row(1).transpose().array()*Opacity.col(j).array()).transpose().matrix();
676-
//std::cout << "TOA" << std::endl;
677669
}
678670

679671
if(clden_species_list == "nan"){
680672
Opacity.col(0) = VectorXd::Ones(stellar_input.row(0).size());
681673
Stellar_activity.col(j) = (stellar_input.row(1).transpose());}
682-
// std::cout << "Setting TOA opacity" << std::endl;
683674
}
684675

685676
for(int rx = 0; rx < PhotoRxn; rx++){
686677
double j_rate = QPhotoChemRate(stellar_input.row(0),d_wavelength, photo_cross_data.col(rx), qyield_data.col(rx), Stellar_activity.col(j));
687678
gas_kin2->setMultiplier(RxnIndex(rx), j_rate);
688679
auto& rxnObj = *(gas_kin->reaction(RxnIndex(rx)));
689680
std::string rxnEquation = rxnObj.equation();
690-
//std::cout<< j << " " << AtmData(iAlt, j)/1e3 << " " << rxnEquation << " " << j_rate << std::endl;
691681

692682
}
693683

@@ -706,7 +696,6 @@ for(int rx = 0; rx < PhotoRxn; rx++){
706696
TChem.col(j) = (isnan(abs(TChem.col(j).array())) ).select(1e20, TChem.col(j));
707697
TChem.col(j) = (isinf(abs(TChem.col(j).array())) ).select(1e20, TChem.col(j));
708698
TChem.col(j) = (TChem.col(j).array().abs() == 0).select(1e20, TChem.col(j));
709-
// std::cout << TChem.col(j).transpose() << std::endl;
710699
//Solve for diffusion terms
711700
if ((j > 0) && (j < nSize-1)) {
712701

@@ -729,47 +718,37 @@ for(int rx = 0; rx < PhotoRxn; rx++){
729718
VectorXd D_prev = handleCustomMolecularDiffusion(PlanetName, gasThermo, AtmData(iPress,j-1), AtmData(iTemp,j-1), mWt);
730719
D_next = (D_next + D_this)/2;
731720
D_prev = (D_prev + D_this)/2;
732-
//std::cout << AtmData(iAlt,j) << std::endl;
733-
//std::cout << D_prev.transpose() << std::endl;
734-
//std::cout << K_prev << std::endl;
735-
//std::cout << "Averaging diffusion coefficient" << std::endl;
736721
double T_next = (AtmData(iTemp,j) + AtmData(iTemp,j+1))/2;
737722
double T_prev = (AtmData(iTemp,j) + AtmData(iTemp,j-1))/2;
738723
double dz = (AtmData(iAlt,j-1) - AtmData(iAlt,j+1))/2;
739724
VectorXd I = VectorXd::Ones(nsp);
740725
//Dynamic Time Scales
741726
VectorXd Diff = AtmData(iKzz, j)*I + D_this;
742727
TDyn.col(j) = ((dz*dz)/Diff.array()).matrix();
743-
//std::cout << "Identity matrix" << std::endl;
744728
VectorXd d_next = ((mm_next*I - mWt).array()*D_next.array()).matrix()*(1/(2*dz*dz))*(g*1e-3*dz/(6.022E23*1.38E-23*T_next));
745729
VectorXd d_prev = ((mm_prev*I - mWt).array()*D_prev.array()).matrix()*(1/(2*dz*dz))*(g*1e-3*dz/(6.022E23*1.38E-23*T_prev));
746730

747731
double T_this = AtmData(iTemp, j);
748732
VectorXd dTdz_prev = alpha*(T_this - T_prev)*2/((T_this + T_prev)*dz);
749733
VectorXd dTdz_next = alpha*(T_next - T_this)*2/((T_this + T_next)*dz);
750734
d_next = d_next + ((dTdz_next.array()*D_next.array()).matrix()*(1/(2*dz*dz)));
751-
d_prev = d_prev + ((dTdz_next.array()*D_next.array()).matrix()*(1/(2*dz*dz)));
735+
d_prev = d_prev + ((dTdz_prev.array()*D_prev.array()).matrix()*(1/(2*dz*dz)));
752736

753-
//std::cout << "d complete" << std::endl;
754737
VectorXd k_next = (((K_next*I) + D_next)/(dz*dz));
755738
VectorXd k_prev = (((K_prev*I) + D_prev)/(dz*dz));
756-
//std::cout << "k complete" << std::endl;
739+
757740
double N_next = AtmData(iNd, j+1);
758741
double N_this = AtmData(iNd, j);
759742
double N_prev = AtmData(iNd, j-1);
760743
double N_n = (N_this + N_next)/2;
761744
double N_p = (N_this + N_prev)/2;
762-
// std::cout << "Averaging densities complete" << std::endl;
763745

764746
//Contribution of U_i
765747
B = ((k_next*N_n/N_this) + (d_next*N_n/N_this) + (k_prev*N_p/N_this) - (d_prev*N_p/N_this)) ;
766-
// std::cout << "U_i complete" << std::endl;
767748
//Contribution of U_i-1
768749
A = -1*((k_prev*N_p/N_prev) + (d_prev*N_p/N_prev));
769-
// std::cout << "U_i-1 complete" << std::endl;
770750
//Contribution of U_i+1
771751
C = -1*((k_next*N_n/N_next) - (d_next*N_n/N_next));
772-
// std::cout << "U_i+1 complete" << std::endl;
773752
//Check for dn/dt
774753
conv.col(j) = m_wdot - (B.array()*N_this*ChemMoleFrac.col(j).array()).matrix() - (A.array()*N_prev*ChemMoleFrac.col(j-1).array()).matrix() - (C.array()*N_prev*ChemMoleFrac.col(j+1).array()).matrix();
775754

@@ -779,7 +758,6 @@ for(int rx = 0; rx < PhotoRxn; rx++){
779758
//Main diagonal terms - Jacobian from chemistry
780759
BlockMatrix.block(j*nsp, j*nsp, nsp, nsp) = ((mat1) - (m_wjac*dt ));
781760
Un = ((mat1) - (m_wjac*dt ))*N_this*mole_frac;
782-
// std::cout << "main diagonal set" << std::endl;
783761
//Main diagonal - diffusion term
784762
for(int sp = 0; sp < nsp; sp++){
785763
Upresent(j*nsp + sp) = (m_wdot(sp)*dt) + Un(sp);
@@ -789,9 +767,8 @@ for(int rx = 0; rx < PhotoRxn; rx++){
789767
//Off diagonal terms - diffusion
790768
BlockMatrix((j)*nsp + sp, (j-1)*nsp + sp) = A(sp)*dt;
791769
BlockMatrix((j)*nsp + sp, (j+1)*nsp + sp) = C(sp)*dt;
792-
}
793-
// std::cout << "Off diagonal terms set" << std::endl;
794-
}
770+
771+
} }
795772

796773
//Upper boundary condition
797774
if(j == 0){
@@ -847,7 +824,7 @@ for(int rx = 0; rx < PhotoRxn; rx++){
847824
VectorXd dTdz_prev = alpha*(T_this - T_prev)*2/((T_this + T_prev)*dz);
848825
VectorXd dTdz_next = alpha*(T_next - T_this)*2/((T_this + T_next)*dz);
849826
d_next = d_next + ((dTdz_next.array()*D_next.array()).matrix()*(1/(2*dz*dz)));
850-
d_prev = d_prev + ((dTdz_next.array()*D_next.array()).matrix()*(1/(2*dz*dz)));
827+
d_prev = d_prev + ((dTdz_prev.array()*D_prev.array()).matrix()*(1/(2*dz*dz)));
851828

852829

853830
//std::cout << "d complete" << std::endl;
@@ -882,7 +859,7 @@ for(int rx = 0; rx < PhotoRxn; rx++){
882859
Un = ((mat1) - (m_wjac*dt ))*N_this*mole_frac;
883860
//Main diagonal - diffusion term
884861
for(int sp = 0; sp < nsp; sp++){
885-
//Dirichlet conditions only for species for whom boundary conditions have been specified
862+
//Dirichlet conditions only for species for whom boundary conditions have been specifies
886863
std::string speciesName = gas->speciesName(sp);
887864
std::string speciesD = pinput->GetOrAddString("upperboundaryMixRat", speciesName, "nan");
888865
if(speciesD != "nan")
@@ -960,48 +937,37 @@ for(int rx = 0; rx < PhotoRxn; rx++){
960937
VectorXd D_this = handleCustomMolecularDiffusion(PlanetName, gasThermo, AtmData(iPress,j), AtmData(iTemp,j), mWt);
961938
VectorXd D_next = handleCustomMolecularDiffusion(PlanetName, gasThermo, AtmData(iPress,j), AtmData(iTemp,j), mWt);
962939
VectorXd D_prev = handleCustomMolecularDiffusion(PlanetName, gasThermo, AtmData(iPress,j-1), AtmData(iTemp,j-1), mWt);
963-
//std::cout << "Diffusion coefficient complete!" << std::endl;
964940
D_next = (D_next + D_this)/2;
965941
D_prev = (D_prev + D_this)/2;
966-
//std::cout << "Averaging diffusion coefficient" << std::endl;
967942
double T_next = (AtmData(iTemp,j) + AtmData(iTemp,j))/2;
968943
double T_prev = (AtmData(iTemp,j) + AtmData(iTemp,j-1))/2;
969944
double dz = (AtmData(iAlt,j-1) - AtmData(iAlt,j));
970945
VectorXd I = VectorXd::Ones(nsp);
971946
//Dynamic Time Scales
972947
VectorXd Diff = AtmData(iKzz, j)*I + D_this;
973948
TDyn.col(j) = ((dz*dz)/Diff.array()).matrix();
974-
//std::cout << "Identity matrix" << std::endl;
975949
VectorXd d_next = ((mm_next*I - mWt).array()*D_next.array()).matrix()*(1/(2*dz*dz))*(g*1e-3*dz/(6.022E23*1.38E-23*T_next));
976950
VectorXd d_prev = ((mm_prev*I - mWt).array()*D_prev.array()).matrix()*(1/(2*dz*dz))*(g*1e-3*dz/(6.022E23*1.38E-23*T_prev));
977-
//std::cout << (d_next - d_prev).transpose() << std::endl;
978951

979952
double T_this = AtmData(iTemp, j);
980953
VectorXd dTdz_prev = alpha*(T_this - T_prev)*2/((T_this + T_prev)*dz);
981954
VectorXd dTdz_next = alpha*(T_next - T_this)*2/((T_this + T_next)*dz);
982955
d_next = d_next + ((dTdz_next.array()*D_next.array()).matrix()*(1/(2*dz*dz)));
983-
d_prev = d_prev + ((dTdz_next.array()*D_next.array()).matrix()*(1/(2*dz*dz)));
956+
d_prev = d_prev + ((dTdz_prev.array()*D_prev.array()).matrix()*(1/(2*dz*dz)));
984957

985-
986-
//std::cout << "d complete" << std::endl;
987958
VectorXd k_next = (((K_next*I) + D_next)/(dz*dz));
988959
VectorXd k_prev = (((K_prev*I) + D_prev)/(dz*dz));
989-
//std::cout << "k complete" << std::endl;
990960
double N_next = AtmData(iNd, j);
991961
double N_this = AtmData(iNd, j);
992962
double N_prev = AtmData(iNd, j-1);
993963
double N_n = (N_this + N_next)/2;
994964
double N_p = (N_this + N_prev)/2;
995-
// std::cout << "Averaging densities complete" << std::endl;
996965
//Contribution of U_i
997966
B = ((k_next*N_n/N_this) + (d_next*N_n/N_this) + (k_prev*N_p/N_this) - (d_prev*N_p/N_this)) ;
998-
// std::cout << "U_i complete" << std::endl;
999967
//Contribution of U_i-1
1000968
A = -1*((k_prev*N_p/N_prev) + (d_prev*N_p/N_prev));
1001-
// std::cout << "U_i-1 complete" << std::endl;
1002969
//Contribution of U_i+1
1003970
C = -1*((k_next*N_n/N_next) - (d_next*N_n/N_next));
1004-
// std::cout << "U_i+1 complete" << std::endl;
1005971

1006972
//Check for dn/dt
1007973
conv.col(j) = m_wdot - (B.array()*N_this*mole_frac.array()).matrix() - (A.array()*N_prev*ChemMoleFrac.col(j-1).array()).matrix() - (C.array()*N_prev*mole_frac.array()).matrix();
@@ -1056,25 +1022,28 @@ for(int rx = 0; rx < PhotoRxn; rx++){
10561022
}
10571023

10581024
//Inverting the matrix
1059-
Eigen::SparseMatrix<double> Am = BlockMatrix.sparseView();
1060-
Eigen::BiCGSTAB<Eigen::SparseMatrix<double>, Eigen::DiagonalPreconditioner<double>> cg;
1061-
cg.setMaxIterations(1E8);
1062-
cg.compute(Am);
1063-
Ufuture = cg.solve(Upresent);
1025+
1026+
// Eigen::SparseMatrix<double> Am = BlockMatrix.sparseView();
1027+
// Eigen::BiCGSTAB<Eigen::SparseMatrix<double>, Eigen::DiagonalPreconditioner<double>> cg;
1028+
// cg.setMaxIterations(1E8);
1029+
// cg.compute(Am);
1030+
// Ufuture = cg.solve(Upresent);
10641031

10651032
// Ufuture = BlockMatrix.inverse()*Upresent;
10661033

1034+
Ufuture = ThomasSolver(BlockMatrix, Upresent,nsp,nSize);
1035+
10671036
//Checking for convergence
10681037
double conv_factor = conv.maxCoeff();
10691038
double tau = 1.1;
1070-
std::cout << "#iterations: " << cg.iterations() << std::endl;
1071-
std::cout << "estimated error: " << cg.error() << std::endl;
1039+
// std::cout << "#iterations: " << cg.iterations() << std::endl;
1040+
// std::cout << "estimated error: " << cg.error() << std::endl;
10721041
std::cout << "convergence: " << conv_factor << std::endl;
10731042
std::cout << "Chemical time scale: " << TChem.minCoeff() << std::endl;
10741043
std::cout << "Dynamic time scale: " << TDyn.minCoeff() << std::endl;
10751044
//Condition number of matrix
1076-
double cd_num = (BlockMatrix.squaredNorm())*(BlockMatrix.inverse().squaredNorm());
1077-
std::cout << "Condition number (L^2): " << cd_num << std::endl;
1045+
// double cd_num = (BlockMatrix.squaredNorm())*(BlockMatrix.inverse().squaredNorm());
1046+
// std::cout << "Condition number (L^2): " << cd_num << std::endl;
10781047

10791048
if(counter >= 1){
10801049
if(conv_factor/conv_old >= 1){

0 commit comments

Comments
 (0)