#include #include #include #include #define xm 3 /* number of elements in horizontal direction */ #define ym 4 /* number of elememts in vertical direction */ #define sstol 1.e-3 #define method 2 /* 1 = Bi-CGSTAB, 2 = GMRES */ #define restart 20 /* for GMRES */ #define pr 0 /* 1 = screen output of each iteration, 0 otherwise */ /******************************************************************************** * differential equation: * * Q(x,y) du/dt + A(x,y,t) d2u/dx2 + B(x,y,t) d2u/dxdy + C(x,y,t) d2u/dy2 * * + D(x,y,t) du/dx + E(x,y,t) du/dy + F(x,y,t) u = H(x,y,t) * ********************************************************************************/ #define A(x,y,t) (1) #define B(x,y,t) (1) #define C(x,y,t) (1) #define D(x,y,t) (1) #define E(x,y,t) (1) #define F(x,y,t) (1) #define Q(x,y) (-0.01*x*y) /************************************************* * exact solution (use for initial conditions) * *************************************************/ #define S(x,y,t) ((x+1) * (x+1) * (y+1) * (y+1) * (1+exp(-t))) #define S_x(x,y,t) (2 * (x+1) * (y+1) * (y+1) * (1+exp(-t))) #define S_y(x,y,t) (2 * (x+1) * (x+1) * (y+1) * (1+exp(-t))) #define S_xy(x,y,t) (4 * (x+1) * (y+1) * (1+exp(-t))) #define S_xx(x,y,t) (2 * (y+1) * (y+1) * (1+exp(-t))) #define S_yy(x,y,t) (2 * (x+1) * (x+1) * (1+exp(-t))) #define S_t(x,y,t) ((x+1) * (x+1) * (y+1) * (y+1) * (-exp(-t))) /********************** * forcing function * **********************/ #define H(x,y,t) (A(x,y,t)*S_xx(x,y,t) + B(x,y,t)*S_xy(x,y,t) + C(x,y,t)*S_yy(x,y,t) + D(x,y,t)*S_x(x,y,t) + E(x,y,t)*S_y(x,y,t) + F(x,y,t)*S(x,y,t) + Q(x,y)*S_t(x,y,t)) /**************************************** * Dirichlet (1) or Neumann (0) BC's * ****************************************/ #define BCTypeS 1 #define BCTypeN 0 #define BCTypeW 1 #define BCTypeE 0 /************************ * boundary condtions * ************************/ #define BCValuS(x,t) (S(x,2,t)) #define BCValudS(x,t) (S_x(x,2,t)) #define BCValuN(x,t) (S_y(x,3,t)) #define BCValudN(x,t) (S_xy(x,3,t)) #define BCValuW(y,t) (S(1,y,t)) #define BCValudW(y,t) (S_y(1,y,t)) #define BCValuE(y,t) (S_x(2,y,t)) #define BCValudE(y,t) (S_xy(2,y,t)) #define hermite_fL(n,h) (0.5 * (1+2*n)*(1+2*n)*(1-n)) #define hermite_fR(n,h) (0.5 * (1-2*n)*(1-2*n)*(1+n)) #define hermite_gL(n,h) (h/8. * (2*n+1)*(2*n+1)*(2*n-1)) #define hermite_gR(n,h) (h/8. * (2*n-1)*(2*n-1)*(2*n+1)) #define hermite_dfL(n,h) (3./(2.*h) * (1+2*n) * (1-2*n)) #define hermite_dfR(n,h) (-hermite_dfL(n,h)) #define hermite_dgL(n,h) (0.25 * (2*n+1) * (6*n-1)) #define hermite_dgR(n,h) (0.25 * (2*n-1) * (6*n+1)) #define hermite_d2fL(n,h) (-12.*n/(h*h)) #define hermite_d2fR(n,h) (-hermite_d2fL(n,h)) #define hermite_d2gL(n,h) ((6*n+1) / h) #define hermite_d2gR(n,h) ((6*n-1) / h) #define L_fLx_fLy(x,y,nx,ny,hx,hy,t) (A(x,y,t)*hermite_d2fL(nx,hx)*hermite_fL(ny,hy) + B(x,y,t)*hermite_dfL(nx,hx)*hermite_dfL(ny,hy) + C(x,y,t)*hermite_fL(nx,hx)*hermite_d2fL(ny,hy) + D(x,y,t)*hermite_dfL(nx,hx)*hermite_fL(ny,hy) + E(x,y,t)*hermite_fL(nx,hx)*hermite_dfL(ny,hy) + F(x,y,t)*hermite_fL(nx,hx)*hermite_fL(ny,hy)) #define L_fLx_fRy(x,y,nx,ny,hx,hy,t) (A(x,y,t)*hermite_d2fL(nx,hx)*hermite_fR(ny,hy) + B(x,y,t)*hermite_dfL(nx,hx)*hermite_dfR(ny,hy) + C(x,y,t)*hermite_fL(nx,hx)*hermite_d2fR(ny,hy) + D(x,y,t)*hermite_dfL(nx,hx)*hermite_fR(ny,hy) + E(x,y,t)*hermite_fL(nx,hx)*hermite_dfR(ny,hy) + F(x,y,t)*hermite_fL(nx,hx)*hermite_fR(ny,hy)) #define L_fLx_gLy(x,y,nx,ny,hx,hy,t) (A(x,y,t)*hermite_d2fL(nx,hx)*hermite_gL(ny,hy) + B(x,y,t)*hermite_dfL(nx,hx)*hermite_dgL(ny,hy) + C(x,y,t)*hermite_fL(nx,hx)*hermite_d2gL(ny,hy) + D(x,y,t)*hermite_dfL(nx,hx)*hermite_gL(ny,hy) + E(x,y,t)*hermite_fL(nx,hx)*hermite_dgL(ny,hy) + F(x,y,t)*hermite_fL(nx,hx)*hermite_gL(ny,hy)) #define L_fLx_gRy(x,y,nx,ny,hx,hy,t) (A(x,y,t)*hermite_d2fL(nx,hx)*hermite_gR(ny,hy) + B(x,y,t)*hermite_dfL(nx,hx)*hermite_dgR(ny,hy) + C(x,y,t)*hermite_fL(nx,hx)*hermite_d2gR(ny,hy) + D(x,y,t)*hermite_dfL(nx,hx)*hermite_gR(ny,hy) + E(x,y,t)*hermite_fL(nx,hx)*hermite_dgR(ny,hy) + F(x,y,t)*hermite_fL(nx,hx)*hermite_gR(ny,hy)) #define L_fRx_fLy(x,y,nx,ny,hx,hy,t) (A(x,y,t)*hermite_d2fR(nx,hx)*hermite_fL(ny,hy) + B(x,y,t)*hermite_dfR(nx,hx)*hermite_dfL(ny,hy) + C(x,y,t)*hermite_fR(nx,hx)*hermite_d2fL(ny,hy) + D(x,y,t)*hermite_dfR(nx,hx)*hermite_fL(ny,hy) + E(x,y,t)*hermite_fR(nx,hx)*hermite_dfL(ny,hy) + F(x,y,t)*hermite_fR(nx,hx)*hermite_fL(ny,hy)) #define L_fRx_fRy(x,y,nx,ny,hx,hy,t) (A(x,y,t)*hermite_d2fR(nx,hx)*hermite_fR(ny,hy) + B(x,y,t)*hermite_dfR(nx,hx)*hermite_dfR(ny,hy) + C(x,y,t)*hermite_fR(nx,hx)*hermite_d2fR(ny,hy) + D(x,y,t)*hermite_dfR(nx,hx)*hermite_fR(ny,hy) + E(x,y,t)*hermite_fR(nx,hx)*hermite_dfR(ny,hy) + F(x,y,t)*hermite_fR(nx,hx)*hermite_fR(ny,hy)) #define L_fRx_gLy(x,y,nx,ny,hx,hy,t) (A(x,y,t)*hermite_d2fR(nx,hx)*hermite_gL(ny,hy) + B(x,y,t)*hermite_dfR(nx,hx)*hermite_dgL(ny,hy) + C(x,y,t)*hermite_fR(nx,hx)*hermite_d2gL(ny,hy) + D(x,y,t)*hermite_dfR(nx,hx)*hermite_gL(ny,hy) + E(x,y,t)*hermite_fR(nx,hx)*hermite_dgL(ny,hy) + F(x,y,t)*hermite_fR(nx,hx)*hermite_gL(ny,hy)) #define L_fRx_gRy(x,y,nx,ny,hx,hy,t) (A(x,y,t)*hermite_d2fR(nx,hx)*hermite_gR(ny,hy) + B(x,y,t)*hermite_dfR(nx,hx)*hermite_dgR(ny,hy) + C(x,y,t)*hermite_fR(nx,hx)*hermite_d2gR(ny,hy) + D(x,y,t)*hermite_dfR(nx,hx)*hermite_gR(ny,hy) + E(x,y,t)*hermite_fR(nx,hx)*hermite_dgR(ny,hy) + F(x,y,t)*hermite_fR(nx,hx)*hermite_gR(ny,hy)) #define L_gLx_fLy(x,y,nx,ny,hx,hy,t) (A(x,y,t)*hermite_d2gL(nx,hx)*hermite_fL(ny,hy) + B(x,y,t)*hermite_dgL(nx,hx)*hermite_dfL(ny,hy) + C(x,y,t)*hermite_gL(nx,hx)*hermite_d2fL(ny,hy) + D(x,y,t)*hermite_dgL(nx,hx)*hermite_fL(ny,hy) + E(x,y,t)*hermite_gL(nx,hx)*hermite_dfL(ny,hy) + F(x,y,t)*hermite_gL(nx,hx)*hermite_fL(ny,hy)) #define L_gLx_fRy(x,y,nx,ny,hx,hy,t) (A(x,y,t)*hermite_d2gL(nx,hx)*hermite_fR(ny,hy) + B(x,y,t)*hermite_dgL(nx,hx)*hermite_dfR(ny,hy) + C(x,y,t)*hermite_gL(nx,hx)*hermite_d2fR(ny,hy) + D(x,y,t)*hermite_dgL(nx,hx)*hermite_fR(ny,hy) + E(x,y,t)*hermite_gL(nx,hx)*hermite_dfR(ny,hy) + F(x,y,t)*hermite_gL(nx,hx)*hermite_fR(ny,hy)) #define L_gLx_gLy(x,y,nx,ny,hx,hy,t) (A(x,y,t)*hermite_d2gL(nx,hx)*hermite_gL(ny,hy) + B(x,y,t)*hermite_dgL(nx,hx)*hermite_dgL(ny,hy) + C(x,y,t)*hermite_gL(nx,hx)*hermite_d2gL(ny,hy) + D(x,y,t)*hermite_dgL(nx,hx)*hermite_gL(ny,hy) + E(x,y,t)*hermite_gL(nx,hx)*hermite_dgL(ny,hy) + F(x,y,t)*hermite_gL(nx,hx)*hermite_gL(ny,hy)) #define L_gLx_gRy(x,y,nx,ny,hx,hy,t) (A(x,y,t)*hermite_d2gL(nx,hx)*hermite_gR(ny,hy) + B(x,y,t)*hermite_dgL(nx,hx)*hermite_dgR(ny,hy) + C(x,y,t)*hermite_gL(nx,hx)*hermite_d2gR(ny,hy) + D(x,y,t)*hermite_dgL(nx,hx)*hermite_gR(ny,hy) + E(x,y,t)*hermite_gL(nx,hx)*hermite_dgR(ny,hy) + F(x,y,t)*hermite_gL(nx,hx)*hermite_gR(ny,hy)) #define L_gRx_fLy(x,y,nx,ny,hx,hy,t) (A(x,y,t)*hermite_d2gR(nx,hx)*hermite_fL(ny,hy) + B(x,y,t)*hermite_dgR(nx,hx)*hermite_dfL(ny,hy) + C(x,y,t)*hermite_gR(nx,hx)*hermite_d2fL(ny,hy) + D(x,y,t)*hermite_dgR(nx,hx)*hermite_fL(ny,hy) + E(x,y,t)*hermite_gR(nx,hx)*hermite_dfL(ny,hy) + F(x,y,t)*hermite_gR(nx,hx)*hermite_fL(ny,hy)) #define L_gRx_fRy(x,y,nx,ny,hx,hy,t) (A(x,y,t)*hermite_d2gR(nx,hx)*hermite_fR(ny,hy) + B(x,y,t)*hermite_dgR(nx,hx)*hermite_dfR(ny,hy) + C(x,y,t)*hermite_gR(nx,hx)*hermite_d2fR(ny,hy) + D(x,y,t)*hermite_dgR(nx,hx)*hermite_fR(ny,hy) + E(x,y,t)*hermite_gR(nx,hx)*hermite_dfR(ny,hy) + F(x,y,t)*hermite_gR(nx,hx)*hermite_fR(ny,hy)) #define L_gRx_gLy(x,y,nx,ny,hx,hy,t) (A(x,y,t)*hermite_d2gR(nx,hx)*hermite_gL(ny,hy) + B(x,y,t)*hermite_dgR(nx,hx)*hermite_dgL(ny,hy) + C(x,y,t)*hermite_gR(nx,hx)*hermite_d2gL(ny,hy) + D(x,y,t)*hermite_dgR(nx,hx)*hermite_gL(ny,hy) + E(x,y,t)*hermite_gR(nx,hx)*hermite_dgL(ny,hy) + F(x,y,t)*hermite_gR(nx,hx)*hermite_gL(ny,hy)) #define L_gRx_gRy(x,y,nx,ny,hx,hy,t) (A(x,y,t)*hermite_d2gR(nx,hx)*hermite_gR(ny,hy) + B(x,y,t)*hermite_dgR(nx,hx)*hermite_dgR(ny,hy) + C(x,y,t)*hermite_gR(nx,hx)*hermite_d2gR(ny,hy) + D(x,y,t)*hermite_dgR(nx,hx)*hermite_gR(ny,hy) + E(x,y,t)*hermite_gR(nx,hx)*hermite_dgR(ny,hy) + F(x,y,t)*hermite_gR(nx,hx)*hermite_gR(ny,hy)) #define J_fLx_fLy(x,y,nx,ny,hx,hy) (Q(x,y) * hermite_fL(nx,hx) * hermite_fL(ny,hy)) #define J_fLx_fRy(x,y,nx,ny,hx,hy) (Q(x,y) * hermite_fL(nx,hx) * hermite_fR(ny,hy)) #define J_fLx_gLy(x,y,nx,ny,hx,hy) (Q(x,y) * hermite_fL(nx,hx) * hermite_gL(ny,hy)) #define J_fLx_gRy(x,y,nx,ny,hx,hy) (Q(x,y) * hermite_fL(nx,hx) * hermite_gR(ny,hy)) #define J_fRx_fLy(x,y,nx,ny,hx,hy) (Q(x,y) * hermite_fR(nx,hx) * hermite_fL(ny,hy)) #define J_fRx_fRy(x,y,nx,ny,hx,hy) (Q(x,y) * hermite_fR(nx,hx) * hermite_fR(ny,hy)) #define J_fRx_gLy(x,y,nx,ny,hx,hy) (Q(x,y) * hermite_fR(nx,hx) * hermite_gL(ny,hy)) #define J_fRx_gRy(x,y,nx,ny,hx,hy) (Q(x,y) * hermite_fR(nx,hx) * hermite_gR(ny,hy)) #define J_gLx_fLy(x,y,nx,ny,hx,hy) (Q(x,y) * hermite_gL(nx,hx) * hermite_fL(ny,hy)) #define J_gLx_fRy(x,y,nx,ny,hx,hy) (Q(x,y) * hermite_gL(nx,hx) * hermite_fR(ny,hy)) #define J_gLx_gLy(x,y,nx,ny,hx,hy) (Q(x,y) * hermite_gL(nx,hx) * hermite_gL(ny,hy)) #define J_gLx_gRy(x,y,nx,ny,hx,hy) (Q(x,y) * hermite_gL(nx,hx) * hermite_gR(ny,hy)) #define J_gRx_fLy(x,y,nx,ny,hx,hy) (Q(x,y) * hermite_gR(nx,hx) * hermite_fL(ny,hy)) #define J_gRx_fRy(x,y,nx,ny,hx,hy) (Q(x,y) * hermite_gR(nx,hx) * hermite_fR(ny,hy)) #define J_gRx_gLy(x,y,nx,ny,hx,hy) (Q(x,y) * hermite_gR(nx,hx) * hermite_gL(ny,hy)) #define J_gRx_gRy(x,y,nx,ny,hx,hy) (Q(x,y) * hermite_gR(nx,hx) * hermite_gR(ny,hy)) #define maxiter 200 #define xc (2 * xm) /* number of collocation points on each row */ #define yc (2 * ym) /* number of collocation points on each column */ #define TRUE 1 #define FALSE 0 #define restol 1.e-8 #define pi (4*atan(1.)) #include "struct2d.h" #include "fun2d_p1.h" #include "fun2d_p2.h" #include "fun2d_p3.h" main() { struct collmat G, T, M; struct collvec U, b, bcombo, t1, cvec, cbar; struct corner_element_bcs BCMSW, BCMSE, BCMNW, BCMNE, BCTSW, BCTSE, BCTNW, BCTNE; struct EW_element_bcs BCMW, BCME, BCTW, BCTE; struct NS_element_bcs BCMN, BCMS, BCTN, BCTS; int unkSW, unkSE, unkNW, unkNE, i, j, k, iter, tsteps, ti; double *x, *y, *Wv, *dWv, *Ev, *dEv, *Sv, *dSv, *Nv, *dNv, nm, **xi, temp, midx, midy, *hx, *hy, **ccx, **ccy, **cx, **cy, t, init_t, dt, theta, tau, taubar, ***SS, **V, **oldV; /****************** * define nodes * ******************/ x = dvector(xm+1); y = dvector(ym+1); x[0] = 1.0; x[1] = 1.3; x[2] = 1.8; x[3] = 2.0; y[0] = 2.0; y[1] = 2.2; y[2] = 2.5; y[3] = 2.9; y[4] = 3.0; /********************** * define xi values * **********************/ xi = dmatrix(xm, ym); for (i = 0; i < xm; i++) for (j = 0; j < ym; j++) xi[i][j] = 1. / sqrt(12); /****************** * time control * ******************/ init_t = 0; dt = 0.1; tsteps = 200; theta = 0.5; tau = theta * dt; taubar = (1-theta) * dt; /********************************************************** * define finite element lengths and collocation points * **********************************************************/ hx = dvector(xm); hy = dvector(ym); cx = dmatrix(xc, yc); cy = dmatrix(xc, yc); ccx = dmatrix(xc, yc); ccy = dmatrix(xc, yc); for (i = 0; i < xm; i++) hx[i] = x[i+1] - x[i]; for (j = 0; j < ym; j++) hy[j] = y[j+1] - y[j]; for (i = 0; i < xm; i++) { midx = 0.5 * (x[i] + x[i+1]); for (j = 0; j < ym; j++) { midy = 0.5 * (y[j] + y[j+1]); cx [2*i] [2*j] = - xi[i][j]; cy [2*i] [2*j] = - xi[i][j]; ccx[2*i] [2*j] = midx - xi[i][j] * hx[i]; ccy[2*i] [2*j] = midy - xi[i][j] * hy[j]; cx [2*i] [2*j+1] = - xi[i][j]; cy [2*i] [2*j+1] = xi[i][j]; ccx[2*i] [2*j+1] = midx - xi[i][j] * hx[i]; ccy[2*i] [2*j+1] = midy + xi[i][j] * hy[j]; cx [2*i+1][2*j] = xi[i][j]; cy [2*i+1][2*j] = - xi[i][j]; ccx[2*i+1][2*j] = midx + xi[i][j] * hx[i]; ccy[2*i+1][2*j] = midy - xi[i][j] * hy[j]; cx [2*i+1][2*j+1] = xi[i][j]; cy [2*i+1][2*j+1] = xi[i][j]; ccx[2*i+1][2*j+1] = midx + xi[i][j] * hx[i]; ccy[2*i+1][2*j+1] = midy + xi[i][j] * hy[j]; } } /******************************* * determine corner unknowns * *******************************/ BCind(&unkNE, &unkNW, &unkSE, &unkSW); /************************ * initial conditions * ************************/ t = init_t; load_collvec_ic(&U, x, y, unkSW, unkSE, unkNW, unkNE, t); /************************************** * time-independent and initial BCs * **************************************/ Wv = dvector(ym+1); Ev = dvector(ym+1); dWv = dvector(ym+1); dEv = dvector(ym+1); Sv = dvector(xm+1); Nv = dvector(xm+1); dSv = dvector(xm+1); dNv = dvector(xm+1); load_BC_J(&BCTSW, &BCTSE, &BCTNW, &BCTNE, &BCTW, &BCTE, &BCTN, &BCTS, hx, hy, cx, cy, ccx, ccy); load_BC_L(&BCMSW, &BCMSE, &BCMNW, &BCMNE, &BCMW, &BCME, &BCMN, &BCMS, hx, hy, cx, cy, ccx, ccy, t); bdy_vals(Wv, dWv, Ev, dEv, Sv, dSv, Nv, dNv, x, y, t); /****************************** * matrices T and initial M * ******************************/ make_T(&T, hx, hy, unkSW, unkSE, unkNW, unkNE, ccx, ccy, cx, cy); make_M(&M, hx, hy, unkSW, unkSE, unkNW, unkNE, ccx, ccy, cx, cy, t); /********************************************** * initial RHS vector at collocation points * **********************************************/ for (i = 0; i < xm; i++) { for (j = 0; j < 2; j++) { b.F.el[i].el[j] = H(ccx[2*i+j][0], ccy[2*i+j][0], t); b.L.el[i].el[j] = H(ccx[2*i+j][2*ym-1], ccy[2*i+j][2*ym-1], t); } for (j = 0; j <= ym-2; j++) { b.M[j].el[i].el[0].el[0] = H(ccx[2*i] [2*j+1], ccy[2*i] [2*j+1] ,t); b.M[j].el[i].el[0].el[1] = H(ccx[2*i] [2*j+2], ccy[2*i] [2*j+2] ,t); b.M[j].el[i].el[1].el[0] = H(ccx[2*i+1][2*j+1], ccy[2*i+1][2*j+1] ,t); b.M[j].el[i].el[1].el[1] = H(ccx[2*i+1][2*j+2], ccy[2*i+1][2*j+2] ,t); } } /****************** * time advance * ******************/ V = dmatrix(ym+1, xm+1); oldV = dmatrix(ym+1, xm+1); SS = d3matrix(4, ym+1, xm+1); for (ti = 1; ti <= tsteps; ti++) { /**************************************** * matrix * solution at old time step * ****************************************/ saxpy_collmat(&G, -taubar, &M, &T); mult_collvec_collmat(&t1, &G, &U); /*************************************** * forcing function at old time step * ***************************************/ scmult_collvec(&cbar, &b, taubar); sub_BC(&cbar, &BCTSW, &BCTSE, &BCTNW, &BCTNE, &BCTW, &BCTE, &BCTN, &BCTS, &BCMSW, &BCMSE, &BCMNW, &BCMNE, &BCMW, &BCME, &BCMN, &BCMS, Wv, Ev, Nv, Sv, dWv, dEv, dNv, dSv, unkSW, unkSE, unkNW, unkNE, taubar, 1.); /******************* * new time step * *******************/ t += dt; /******************* * make matrix G * *******************/ make_M(&M, hx, hy, unkSW, unkSE, unkNW, unkNE, ccx, ccy, cx, cy, t); saxpy_collmat(&G, tau, &M, &T); /************************** * BCs at new time step * **************************/ load_BC_L(&BCMSW, &BCMSE, &BCMNW, &BCMNE, &BCMW, &BCME, &BCMN, &BCMS, hx, hy, cx, cy, ccx, ccy, t); bdy_vals(Wv, dWv, Ev, dEv, Sv, dSv, Nv, dNv, x, y, t); /*********************************************** * evaluate RHS vector at collocation points * ***********************************************/ for (i = 0; i < xm; i++) { for (j = 0; j < 2; j++) { b.F.el[i].el[j] = H(ccx[2*i+j][0], ccy[2*i+j][0], t); b.L.el[i].el[j] = H(ccx[2*i+j][2*ym-1], ccy[2*i+j][2*ym-1], t); } for (j = 0; j <= ym-2; j++) { b.M[j].el[i].el[0].el[0] = H(ccx[2*i] [2*j+1], ccy[2*i] [2*j+1] ,t); b.M[j].el[i].el[0].el[1] = H(ccx[2*i] [2*j+2], ccy[2*i] [2*j+2] ,t); b.M[j].el[i].el[1].el[0] = H(ccx[2*i+1][2*j+1], ccy[2*i+1][2*j+1] ,t); b.M[j].el[i].el[1].el[1] = H(ccx[2*i+1][2*j+2], ccy[2*i+1][2*j+2] ,t); } } scmult_collvec(&cvec, &b, tau); /********************************** * subtract BCs from RHS vector * **********************************/ sub_BC(&cvec, &BCTSW, &BCTSE, &BCTNW, &BCTNE, &BCTW, &BCTE, &BCTN, &BCTS, &BCMSW, &BCMSE, &BCMNW, &BCMNE, &BCMW, &BCME, &BCMN, &BCMS, Wv, Ev, Nv, Sv, dWv, dEv, dNv, dSv, unkSW, unkSE, unkNW, unkNE, tau, -1.); /*************************** * form final RHS vector * ***************************/ add_collvec(&bcombo, &cvec, &cbar); add_collvec(&t1, &t1, &bcombo); /************************************** * compute solution vector via PBCG * **************************************/ iter = 0; if (method == 1) pbcgstab(&U, &G, &t1, &iter, &nm); else gmres(&U, &G, &t1, &iter, &nm); /*********************** * results to screen * ***********************/ soln_gen(SS, &U, Wv, dWv, Ev, dEv, Sv, dSv, Nv, dNv, unkSW, unkSE, unkNW, unkNE); copymat(oldV, V); copymat(V, SS[0]); blank(1); printf("time step = %d\ntime = %lf\ndt = %lf\n", ti, t, dt); for (i = 0; i <= xm; i++) { blank(1); for (j = ym; j >= 0; j--) printf("%8.3lf %8.3lf %8.3lf %8.3lf %8.3lf %8.3lf\n", x[i], y[ym-j], SS[0][j][i], SS[1][j][i], SS[2][j][i], SS[3][j][i]); } temp = compmat(oldV, V); printf("max = %lf\n", temp); if (temp < sstol) { printf("steady state reached\n"); ti = tsteps; } if (iter == maxiter) { printf("\nmaximum number of %d iterations reached\nExiting program...\n\n", iter); exit(0); } blank(1); } } /* end of main */