| c--------------------------------------------------------- |
| subroutine convect(ifmortar) |
| c--------------------------------------------------------- |
| c Advance the convection term using 4th order RK |
| c 1.ta1 is solution from last time step |
| c 2.the heat source is considered part of d/dx |
| c 3.trhs is right hand side for the diffusion equation |
| c 4.tmor is solution on mortar points, which will be used |
| c as the initial guess when advancing the diffusion term |
| c--------------------------------------------------------- |
| |
| include 'header.h' |
| |
| double precision alpha2, tempa(lx1,lx1,lx1), |
| & rdtime, pidivalpha, sixth, |
| & dtx1, dtx2, dtx3, src, rk1(lx1,lx1,lx1), rk2(lx1,lx1,lx1), |
| & rk3(lx1,lx1,lx1), rk4(lx1,lx1,lx1), temp(lx1,lx1,lx1), |
| & subtime(3), xx0(3), yy0(3), zz0(3), dtime2, r2, sum, |
| & xloc(lx1), yloc(lx1), zloc(lx1) |
| integer k,iel,i,j,iside,isize, substep, ip |
| logical ifmortar |
| parameter (sixth=1.d0/6.d0) |
| |
| if (timeron) call timer_start(t_convect) |
| pidivalpha = dacos(-1.d0)/alpha |
| alpha2 = alpha*alpha |
| dtime2 = dtime/2.d0 |
| rdtime = 1.d0/dtime |
| subtime(1) = time |
| subtime(2) = time+dtime2 |
| subtime(3) = time+dtime |
| do substep = 1, 3 |
| xx0(substep) = x00+velx*subtime(substep) |
| yy0(substep) = y00+vely*subtime(substep) |
| zz0(substep) = z00+velz*subtime(substep) |
| end do |
| |
| c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(rk4,rk3,rk2,temp,rk1,dtx3, |
| c$OMP$ dtx2,dtx1,iside,ip,sum,src,r2,i,j,k,isize,iel,tempa, |
| c$OMP$ xloc,yloc,zloc) |
| |
| do iel = 1, nelt |
| isize=size_e(iel) |
| c.......xloc(i) is the location of i'th collocation in x direction in an element. |
| c yloc(i) is the location of j'th collocation in y direction in an element. |
| c zloc(i) is the location of k'th collocation in z direction in an element. |
| do i = 1, lx1 |
| xloc(i) = xfrac(i)*(xc(2,iel)-xc(1,iel))+xc(1,iel) |
| end do |
| do j = 1, lx1 |
| yloc(j) = xfrac(j)*(yc(4,iel)-yc(1,iel))+yc(1,iel) |
| end do |
| do k = 1, lx1 |
| zloc(k) = xfrac(k)*(zc(5,iel)-zc(1,iel))+zc(1,iel) |
| end do |
| |
| do k = 1, lx1 |
| do j = 1, lx1 |
| do i = 1, lx1 |
| r2 = (xloc(i)-xx0(1))**2+(yloc(j)-yy0(1))**2+ |
| & (zloc(k)-zz0(1))**2 |
| if (r2.le.alpha2) then |
| src = dcos(dsqrt(r2)*pidivalpha)+1.d0 |
| else |
| src = 0.d0 |
| endif |
| sum = 0.d0 |
| do ip = 1, lx1 |
| sum = sum + dxm1(i,ip) * ta1(ip,j,k,iel) |
| end do |
| dtx1 = -velx*sum*xrm1_s(i,j,k,isize) |
| sum = 0.d0 |
| do ip = 1, lx1 |
| sum = sum + dxm1(j,ip) * ta1(i,ip,k,iel) |
| end do |
| dtx2=-vely*sum*xrm1_s(i,j,k,isize) |
| sum = 0.d0 |
| do ip = 1, lx1 |
| sum = sum + dxm1(k,ip) * ta1(i,j,ip,iel) |
| end do |
| dtx3=-velz*sum*xrm1_s(i,j,k,isize) |
| |
| rk1(i,j,k)= dtx1 + dtx2 + dtx3 + src |
| temp(i,j,k)=ta1(i,j,k,iel)+dtime2*rk1(i,j,k) |
| |
| end do |
| end do |
| end do |
| |
| do k = 1, lx1 |
| do j = 1, lx1 |
| do i = 1, lx1 |
| r2 = (xloc(i)-xx0(2))**2 + (yloc(j)-yy0(2))**2 + |
| & (zloc(k)-zz0(2))**2 |
| if (r2.le.alpha2) then |
| src = dcos(dsqrt(r2)*pidivalpha)+1.d0 |
| else |
| src = 0.d0 |
| endif |
| sum = 0.d0 |
| do ip = 1, lx1 |
| sum = sum + dxm1(i,ip) * temp(ip,j,k) |
| end do |
| dtx1 =-velx*sum*xrm1_s(i,j,k,isize) |
| sum = 0.d0 |
| do ip = 1, lx1 |
| sum = sum + dxm1(j,ip) * temp(i,ip,k) |
| end do |
| dtx2 =-vely*sum*xrm1_s(i,j,k,isize) |
| sum = 0.d0 |
| do ip = 1, lx1 |
| sum = sum + dxm1(k,ip) * temp(i,j,ip) |
| end do |
| dtx3 =-velz*sum*xrm1_s(i,j,k,isize) |
| |
| rk2(i,j,k)= dtx1 + dtx2 + dtx3 + src |
| tempa(i,j,k)=ta1(i,j,k,iel)+dtime2*rk2(i,j,k) |
| end do |
| end do |
| end do |
| |
| do k = 1, lx1 |
| do j = 1, lx1 |
| do i = 1, lx1 |
| r2 = (xloc(i)-xx0(2))**2 + (yloc(j)-yy0(2))**2 + |
| & (zloc(k)-zz0(2))**2 |
| if (r2.le.alpha2) then |
| src = dcos(dsqrt(r2)*pidivalpha)+1.d0 |
| else |
| src = 0.d0 |
| endif |
| sum = 0.d0 |
| do ip = 1, lx1 |
| sum = sum + dxm1(i,ip) * tempa(ip,j,k) |
| end do |
| dtx1 =-velx*sum*xrm1_s(i,j,k,isize) |
| sum = 0.d0 |
| do ip = 1, lx1 |
| sum = sum + dxm1(j,ip) * tempa(i,ip,k) |
| end do |
| dtx2 =-vely*sum*xrm1_s(i,j,k,isize) |
| sum = 0.d0 |
| do ip = 1, lx1 |
| sum = sum + dxm1(k,ip) * tempa(i,j,ip) |
| end do |
| dtx3 =-velz*sum*xrm1_s(i,j,k,isize) |
| |
| rk3(i,j,k)= dtx1 + dtx2 + dtx3 + src |
| temp(i,j,k)=ta1(i,j,k,iel)+dtime*rk3(i,j,k) |
| end do |
| end do |
| end do |
| |
| do k = 1, lx1 |
| do j = 1, lx1 |
| do i = 1, lx1 |
| r2 = (xloc(i)-xx0(3))**2 + (yloc(j)-yy0(3))**2 + |
| & (zloc(k)-zz0(3))**2 |
| if (r2.le.alpha2) then |
| src = dcos(dsqrt(r2)*pidivalpha)+1.d0 |
| else |
| src = 0.d0 |
| endif |
| sum = 0.d0 |
| do ip = 1, lx1 |
| sum = sum + dxm1(i,ip) * temp(ip,j,k) |
| end do |
| dtx1 =-velx*sum*xrm1_s(i,j,k,isize) |
| sum = 0.d0 |
| do ip = 1, lx1 |
| sum = sum + dxm1(j,ip) * temp(i,ip,k) |
| end do |
| dtx2 =-vely*sum*xrm1_s(i,j,k,isize) |
| sum = 0.d0 |
| do ip = 1, lx1 |
| sum = sum + dxm1(k,ip) * temp(i,j,ip) |
| end do |
| dtx3 =-velz*sum*xrm1_s(i,j,k,isize) |
| |
| rk4(i,j,k)= dtx1 + dtx2 + dtx3 + src |
| tempa(i,j,k)=sixth*(rk1(i,j,k)+2.d0* |
| & rk2(i,j,k)+2.d0*rk3(i,j,k)+rk4(i,j,k)) |
| end do |
| end do |
| end do |
| |
| c.......apply boundary condition |
| do iside=1,nsides |
| if(cbc(iside,iel).eq.0)then |
| call facev(tempa,iside,0.0d0) |
| end if |
| end do |
| |
| do k=1,lx1 |
| do j=1,lx1 |
| do i=1,lx1 |
| trhs(i,j,k,iel)=bm1_s(i,j,k,isize)*(ta1(i,j,k,iel)*rdtime+ |
| & tempa(i,j,k)) |
| ta1(i,j,k,iel)=ta1(i,j,k,iel)+tempa(i,j,k)*dtime |
| end do |
| end do |
| end do |
| |
| end do |
| c$OMP END PARALLEL DO |
| |
| c.....get mortar for intial guess for CG |
| |
| if (timeron) call timer_start(t_transfb_c) |
| if(ifmortar)then |
| call transfb_c_2(ta1) |
| else |
| call transfb_c(ta1) |
| end if |
| if (timeron) call timer_stop(t_transfb_c) |
| |
| c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) |
| do i=1,nmor |
| tmort(i)=tmort(i)/mormult(i) |
| end do |
| c$OMP END PARALLEL DO |
| if (timeron) call timer_stop(t_convect) |
| |
| return |
| end |
| |
