/* Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University Copyright Equnior ASA This file is part of the Open Porous Media project (OPM). OPM 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. OPM 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. You should have received a copy of the GNU General Public License along with OPM. If not, see . */ /* Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University Copyright Equnior ASA This file is part of the Open Porous Media project (OPM). OPM 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. OPM 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. You should have received a copy of the GNU General Public License along with OPM. If not, see . */ #include #include extern "C" void ScaLBL_Gradient_Unpack(double weight, double Cqx, double Cqy, double Cqz, int *list, int start, int count, double *recvbuf, double *phi, double *grad, int N){ //.................................................................................... // unpack halo and incorporate into D3Q19 based gradient // Distribution q matche Cqx, Cqy, Cqz //.................................................................................... int n,idx; double value; for (idx=0; idx 0.f){ nA = 1.0; nB = 0.f; } else{ nB = 1.0; nA = 0.f; } Den[idx] = nA; Den[Np+idx] = nB; Aq[idx]=0.3333333333333333*nA; Aq[Np+idx]=0.1111111111111111*nA; Aq[2*Np+idx]=0.1111111111111111*nA; Aq[3*Np+idx]=0.1111111111111111*nA; Aq[4*Np+idx]=0.1111111111111111*nA; Aq[5*Np+idx]=0.1111111111111111*nA; Aq[6*Np+idx]=0.1111111111111111*nA; Bq[idx]=0.3333333333333333*nB; Bq[Np+idx]=0.1111111111111111*nB; Bq[2*Np+idx]=0.1111111111111111*nB; Bq[3*Np+idx]=0.1111111111111111*nB; Bq[4*Np+idx]=0.1111111111111111*nB; Bq[5*Np+idx]=0.1111111111111111*nB; Bq[6*Np+idx]=0.1111111111111111*nB; } } // LBM based on density functional hydrodynamics extern "C" void ScaLBL_D3Q19_AAeven_DFH(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, double *Phi, double *Gradient, double *SolidForce, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int start, int finish, int Np) { double fq; // conserved momemnts double rho,jx,jy,jz; // non-conserved moments double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; double nA,nB; // number density double a1,b1,a2,b2,nAB,delta; double C,nx,ny,nz; //color gradient magnitude and direction double ux,uy,uz; double phi,tau,rho0,rlx_setA,rlx_setB; double force_x,force_y,force_z; const double mrt_V1=0.05263157894736842; const double mrt_V2=0.012531328320802; const double mrt_V3=0.04761904761904762; const double mrt_V4=0.004594820384294068; const double mrt_V5=0.01587301587301587; const double mrt_V6=0.0555555555555555555555555; const double mrt_V7=0.02777777777777778; const double mrt_V8=0.08333333333333333; const double mrt_V9=0.003341687552213868; const double mrt_V10=0.003968253968253968; const double mrt_V11=0.01388888888888889; const double mrt_V12=0.04166666666666666; for (int n=start; n0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta; b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta; Aq[1*Np+n] = a1; Bq[1*Np+n] = b1; Aq[2*Np+n] = a2; Bq[2*Np+n] = b2; //............................................... // q = 2 // Cq = {0,1,0} delta = beta*nA*nB*nAB*0.1111111111111111*ny; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta; b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta; Aq[3*Np+n] = a1; Bq[3*Np+n] = b1; Aq[4*Np+n] = a2; Bq[4*Np+n] = b2; //............................................... // q = 4 // Cq = {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nz; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta; b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta; Aq[5*Np+n] = a1; Bq[5*Np+n] = b1; Aq[6*Np+n] = a2; Bq[6*Np+n] = b2; //............................................... } } extern "C" void ScaLBL_D3Q19_AAodd_DFH(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, double *Phi, double *Gradient, double *SolidForce, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int start, int finish, int Np){ int nread; int nr1,nr2,nr3,nr4,nr5,nr6; int nr7,nr8,nr9,nr10; int nr11,nr12,nr13,nr14; //int nr15,nr16,nr17,nr18; double fq; // conserved momemnts double rho,jx,jy,jz; // non-conserved moments double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; double nA,nB; // number density double a1,b1,a2,b2,nAB,delta; double C,nx,ny,nz; //color gradient magnitude and direction double ux,uy,uz; double phi,tau,rho0,rlx_setA,rlx_setB; double force_x,force_y,force_z; const double mrt_V1=0.05263157894736842; const double mrt_V2=0.012531328320802; const double mrt_V3=0.04761904761904762; const double mrt_V4=0.004594820384294068; const double mrt_V5=0.01587301587301587; const double mrt_V6=0.0555555555555555555555555; const double mrt_V7=0.02777777777777778; const double mrt_V8=0.08333333333333333; const double mrt_V9=0.003341687552213868; const double mrt_V10=0.003968253968253968; const double mrt_V11=0.01388888888888889; const double mrt_V12=0.04166666666666666; for (int n=start; n even part of dist) //fq = dist[nread]; // reading the f2 data into register fq nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist) fq = dist[nr2]; // reading the f2 data into register fq rho += fq; m1 -= 11.0*(fq); m2 -= 4.0*(fq); jx -= fq; m4 += 4.0*(fq); m9 += 2.0*(fq); m10 -= 4.0*(fq); // q=3 //nread = neighborList[n+2*Np]; // neighbor 4 //fq = dist[nread]; nr3 = neighborList[n+2*Np]; // neighbor 4 fq = dist[nr3]; rho += fq; m1 -= 11.0*fq; m2 -= 4.0*fq; jy = fq; m6 = -4.0*fq; m9 -= fq; m10 += 2.0*fq; m11 = fq; m12 = -2.0*fq; // q = 4 //nread = neighborList[n+3*Np]; // neighbor 3 //fq = dist[nread]; nr4 = neighborList[n+3*Np]; // neighbor 3 fq = dist[nr4]; rho+= fq; m1 -= 11.0*fq; m2 -= 4.0*fq; jy -= fq; m6 += 4.0*fq; m9 -= fq; m10 += 2.0*fq; m11 += fq; m12 -= 2.0*fq; // q=5 //nread = neighborList[n+4*Np]; //fq = dist[nread]; nr5 = neighborList[n+4*Np]; fq = dist[nr5]; rho += fq; m1 -= 11.0*fq; m2 -= 4.0*fq; jz = fq; m8 = -4.0*fq; m9 -= fq; m10 += 2.0*fq; m11 -= fq; m12 += 2.0*fq; // q = 6 //nread = neighborList[n+5*Np]; //fq = dist[nread]; nr6 = neighborList[n+5*Np]; fq = dist[nr6]; rho+= fq; m1 -= 11.0*fq; m2 -= 4.0*fq; jz -= fq; m8 += 4.0*fq; m9 -= fq; m10 += 2.0*fq; m11 -= fq; m12 += 2.0*fq; // q=7 //nread = neighborList[n+6*Np]; //fq = dist[nread]; nr7 = neighborList[n+6*Np]; fq = dist[nr7]; rho += fq; m1 += 8.0*fq; m2 += fq; jx += fq; m4 += fq; jy += fq; m6 += fq; m9 += fq; m10 += fq; m11 += fq; m12 += fq; m13 = fq; m16 = fq; m17 = -fq; // q = 8 //nread = neighborList[n+7*Np]; //fq = dist[nread]; nr8 = neighborList[n+7*Np]; fq = dist[nr8]; rho += fq; m1 += 8.0*fq; m2 += fq; jx -= fq; m4 -= fq; jy -= fq; m6 -= fq; m9 += fq; m10 += fq; m11 += fq; m12 += fq; m13 += fq; m16 -= fq; m17 += fq; // q=9 //nread = neighborList[n+8*Np]; //fq = dist[nread]; nr9 = neighborList[n+8*Np]; fq = dist[nr9]; rho += fq; m1 += 8.0*fq; m2 += fq; jx += fq; m4 += fq; jy -= fq; m6 -= fq; m9 += fq; m10 += fq; m11 += fq; m12 += fq; m13 -= fq; m16 += fq; m17 += fq; // q = 10 //nread = neighborList[n+9*Np]; //fq = dist[nread]; nr10 = neighborList[n+9*Np]; fq = dist[nr10]; rho += fq; m1 += 8.0*fq; m2 += fq; jx -= fq; m4 -= fq; jy += fq; m6 += fq; m9 += fq; m10 += fq; m11 += fq; m12 += fq; m13 -= fq; m16 -= fq; m17 -= fq; // q=11 //nread = neighborList[n+10*Np]; //fq = dist[nread]; nr11 = neighborList[n+10*Np]; fq = dist[nr11]; rho += fq; m1 += 8.0*fq; m2 += fq; jx += fq; m4 += fq; jz += fq; m8 += fq; m9 += fq; m10 += fq; m11 -= fq; m12 -= fq; m15 = fq; m16 -= fq; m18 = fq; // q=12 //nread = neighborList[n+11*Np]; //fq = dist[nread]; nr12 = neighborList[n+11*Np]; fq = dist[nr12]; rho += fq; m1 += 8.0*fq; m2 += fq; jx -= fq; m4 -= fq; jz -= fq; m8 -= fq; m9 += fq; m10 += fq; m11 -= fq; m12 -= fq; m15 += fq; m16 += fq; m18 -= fq; // q=13 //nread = neighborList[n+12*Np]; //fq = dist[nread]; nr13 = neighborList[n+12*Np]; fq = dist[nr13]; rho += fq; m1 += 8.0*fq; m2 += fq; jx += fq; m4 += fq; jz -= fq; m8 -= fq; m9 += fq; m10 += fq; m11 -= fq; m12 -= fq; m15 -= fq; m16 -= fq; m18 -= fq; // q=14 //nread = neighborList[n+13*Np]; //fq = dist[nread]; nr14 = neighborList[n+13*Np]; fq = dist[nr14]; rho += fq; m1 += 8.0*fq; m2 += fq; jx -= fq; m4 -= fq; jz += fq; m8 += fq; m9 += fq; m10 += fq; m11 -= fq; m12 -= fq; m15 -= fq; m16 += fq; m18 += fq; // q=15 nread = neighborList[n+14*Np]; fq = dist[nread]; //fq = dist[17*Np+n]; rho += fq; m1 += 8.0*fq; m2 += fq; jy += fq; m6 += fq; jz += fq; m8 += fq; m9 -= 2.0*fq; m10 -= 2.0*fq; m14 = fq; m17 += fq; m18 -= fq; // q=16 nread = neighborList[n+15*Np]; fq = dist[nread]; //fq = dist[8*Np+n]; rho += fq; m1 += 8.0*fq; m2 += fq; jy -= fq; m6 -= fq; jz -= fq; m8 -= fq; m9 -= 2.0*fq; m10 -= 2.0*fq; m14 += fq; m17 -= fq; m18 += fq; // q=17 //fq = dist[18*Np+n]; nread = neighborList[n+16*Np]; fq = dist[nread]; rho += fq; m1 += 8.0*fq; m2 += fq; jy += fq; m6 += fq; jz -= fq; m8 -= fq; m9 -= 2.0*fq; m10 -= 2.0*fq; m14 -= fq; m17 += fq; m18 += fq; // q=18 nread = neighborList[n+17*Np]; fq = dist[nread]; //fq = dist[9*Np+n]; rho += fq; m1 += 8.0*fq; m2 += fq; jy -= fq; m6 -= fq; jz += fq; m8 += fq; m9 -= 2.0*fq; m10 -= 2.0*fq; m14 -= fq; m17 -= fq; m18 -= fq; //........................................................................ //..............carry out relaxation process.............................. //..........Toelke, Fruediger et. al. 2006................................ if (C == 0.0) nx = ny = nz = 0.0; m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) -alpha*C - m1); m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0)- m2); m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4); m6 = m6 + rlx_setB*((-0.6666666666666666*jy)- m6); m8 = m8 + rlx_setB*((-0.6666666666666666*jz)- m8); m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho0) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9); m10 = m10 + rlx_setA*( - m10); m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho0) + 0.5*alpha*C*(ny*ny-nz*nz)- m11); m12 = m12 + rlx_setA*( - m12); m13 = m13 + rlx_setA*( (jx*jy/rho0) + 0.5*alpha*C*nx*ny - m13); m14 = m14 + rlx_setA*( (jy*jz/rho0) + 0.5*alpha*C*ny*nz - m14); m15 = m15 + rlx_setA*( (jx*jz/rho0) + 0.5*alpha*C*nx*nz - m15); m16 = m16 + rlx_setB*( - m16); m17 = m17 + rlx_setB*( - m17); m18 = m18 + rlx_setB*( - m18); // assign force with wetting BC force_x = alpha*(nA-nB)*SolidForce[n] + Fx; force_y = alpha*(nA-nB)*SolidForce[n+Np] + Fy; force_z = alpha*(nA-nB)*SolidForce[n+2*Np] + Fz; //.................inverse transformation...................................................... // q=0 fq = mrt_V1*rho-mrt_V2*m1+mrt_V3*m2; dist[n] = fq; // q = 1 fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10)+0.16666666*force_x; //nread = neighborList[n+Np]; dist[nr2] = fq; // q=2 fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10) - 0.16666666*force_x; //nread = neighborList[n]; dist[nr1] = fq; // q = 3 fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12) + 0.16666666*force_y; //nread = neighborList[n+3*Np]; dist[nr4] = fq; // q = 4 fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12) - 0.16666666*force_y; //nread = neighborList[n+2*Np]; dist[nr3] = fq; // q = 5 fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11) + 0.16666666*force_z; //nread = neighborList[n+5*Np]; dist[nr6] = fq; // q = 6 fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11) - 0.16666666*force_z; //nread = neighborList[n+4*Np]; dist[nr5] = fq; // q = 7 fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)+ mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m16-m17) + 0.08333333333*(force_x+force_y); //nread = neighborList[n+7*Np]; dist[nr8] = fq; // q = 8 fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11 +mrt_V12*m12+0.25*m13+0.125*(m17-m16) - 0.08333333333*(force_x+force_y); //nread = neighborList[n+6*Np]; dist[nr7] = fq; // q = 9 fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)+ mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13+0.125*(m16+m17) + 0.08333333333*(force_x-force_y); //nread = neighborList[n+9*Np]; dist[nr10] = fq; // q = 10 fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)+ mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13-0.125*(m16+m17)- 0.08333333333*(force_x-force_y); //nread = neighborList[n+8*Np]; dist[nr9] = fq; // q = 11 fq = mrt_V1*rho+mrt_V9*m1 +mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8) +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 -mrt_V12*m12+0.25*m15+0.125*(m18-m16) + 0.08333333333*(force_x+force_z); //nread = neighborList[n+11*Np]; dist[nr12] = fq; // q = 12 fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)+ mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m16-m18) - 0.08333333333*(force_x+force_z); //nread = neighborList[n+10*Np]; dist[nr11]= fq; // q = 13 fq = mrt_V1*rho+mrt_V9*m1 +mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8) +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 -mrt_V12*m12-0.25*m15-0.125*(m16+m18) + 0.08333333333*(force_x-force_z); //nread = neighborList[n+13*Np]; dist[nr14] = fq; // q= 14 fq = mrt_V1*rho+mrt_V9*m1 +mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4) +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 -mrt_V12*m12-0.25*m15+0.125*(m16+m18) - 0.08333333333*(force_x-force_z); //nread = neighborList[n+12*Np]; dist[nr13] = fq; // q = 15 fq = mrt_V1*rho+mrt_V9*m1 +mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8) -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18) + 0.08333333333*(force_y+force_z); nread = neighborList[n+15*Np]; dist[nread] = fq; // q = 16 fq = mrt_V1*rho+mrt_V9*m1 +mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8) -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17)- 0.08333333333*(force_y+force_z); nread = neighborList[n+14*Np]; dist[nread] = fq; // q = 17 fq = mrt_V1*rho+mrt_V9*m1 +mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8) -mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18) + 0.08333333333*(force_y-force_z); nread = neighborList[n+17*Np]; dist[nread] = fq; // q = 18 fq = mrt_V1*rho+mrt_V9*m1 +mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6) -mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18) - 0.08333333333*(force_y-force_z); nread = neighborList[n+16*Np]; dist[nread] = fq; // write the velocity ux = (jx + force_x) / rho0; uy = (jy + force_y) / rho0; uz = (jz + force_z) / rho0; //Velocity[n] = ux; //Velocity[Np+n] = uy; //Velocity[2*Np+n] = uz; // Instantiate mass transport distributions // Stationary value - distribution 0 nAB = 1.0/(nA+nB); Aq[n] = 0.3333333333333333*nA; Bq[n] = 0.3333333333333333*nB; //............................................... // q = 0,2,4 // Cq = {1,0,0}, {0,1,0}, {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nx; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta; b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta; // q = 1 //nread = neighborList[n+Np]; Aq[nr2] = a1; Bq[nr2] = b1; // q=2 //nread = neighborList[n]; Aq[nr1] = a2; Bq[nr1] = b2; //............................................... // Cq = {0,1,0} delta = beta*nA*nB*nAB*0.1111111111111111*ny; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta; b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta; // q = 3 //nread = neighborList[n+3*Np]; Aq[nr4] = a1; Bq[nr4] = b1; // q = 4 //nread = neighborList[n+2*Np]; Aq[nr3] = a2; Bq[nr3] = b2; //............................................... // q = 4 // Cq = {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nz; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta; b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta; // q = 5 //nread = neighborList[n+5*Np]; Aq[nr6] = a1; Bq[nr6] = b1; // q = 6 //nread = neighborList[n+4*Np]; Aq[nr5] = a2; Bq[nr5] = b2; //............................................... } } extern "C" void ScaLBL_D3Q7_AAodd_DFH(int *neighborList, double *Aq, double *Bq, double *Den, double *Phi, int start, int finish, int Np) { for (int n=start; n