diff --git a/.travis.yml b/.travis.yml index 18d2ef552..64653dddd 100644 --- a/.travis.yml +++ b/.travis.yml @@ -38,8 +38,8 @@ env: - TEST_CASE=MvCylCvode::test_PnPn_Parallel_Steps1e3 IFMPI=true F77=mpif77 CC=mpicc PPLIST="$PPLIST CVODE" - - TEST_CASE=LinChan_Dir::test_PnPn2_Parallel IFMPI=true F77=mpif77 CC=mpicc - - TEST_CASE=LinChan_Adj::test_PnPn2_Parallel IFMPI=true F77=mpif77 CC=mpicc + - TEST_CASE=LinCav_Dir::test_PnPn2_Parallel IFMPI=true F77=mpif77 CC=mpicc + - TEST_CASE=LinCav_Adj::test_PnPn2_Parallel IFMPI=true F77=mpif77 CC=mpicc before_install: - export ROOT_DIR=`pwd` diff --git a/short_tests/NekTests.py b/short_tests/NekTests.py index d633ad7a7..0137aa0c2 100755 --- a/short_tests/NekTests.py +++ b/short_tests/NekTests.py @@ -939,19 +939,19 @@ def tearDown(self): self.move_logs() #################################################################### -# channel2D; lin_chan_dir.par, lin_chan_adj.par +# dfh_cav; lin_dfh_cav_dir.par, lin_dfh_cav_adj.par #################################################################### -class LinChan_Dir(NekTestCase): - example_subdir = 'channel2D' - case_name = 'lin_chan_dir' +class LinCav_Dir(NekTestCase): + example_subdir = 'dfh_cav' + case_name = 'lin_dfh_cav_dir' def setUp(self): self.size_params = dict( ldim = '2', - lx1 = '10', - lxd = '15', + lx1 = '9', + lxd = '13', lx2 = 'lx1-2', - lelg = '50', + lelg = '500', lpelt = 'lelt', ) @@ -962,11 +962,11 @@ def setUp(self): def test_PnPn2_Serial(self): self.size_params['lx2'] = 'lx1-2' self.config_size() - self.build_nek(usr_file='lin_chan') + self.build_nek(usr_file='lin_dfh_cav') self.run_nek(step_limit=None) omega = self.get_value_from_log('Energy', column=-3, row=-1) - self.assertAlmostEqualDelayed(omega, target_val=-1.2337E-03, delta=2E-06, label='growth rate') + self.assertAlmostEqualDelayed(omega, target_val=-7.57304E-03, delta=1E-06, label='growth rate') self.assertDelayedFailures() @@ -974,28 +974,28 @@ def test_PnPn2_Serial(self): def test_PnPn2_Parallel(self): self.size_params['lx2'] = 'lx1-2' self.config_size() - self.build_nek(usr_file='lin_chan') + self.build_nek(usr_file='lin_dfh_cav') self.run_nek(step_limit=None) omega = self.get_value_from_log('Energy', column=-3, row=-1) - self.assertAlmostEqualDelayed(omega, target_val=-1.2337E-03, delta=2E-06, label='growth rate') + self.assertAlmostEqualDelayed(omega, target_val=-7.57304E-03, delta=1E-06, label='growth rate') self.assertDelayedFailures() def tearDown(self): self.move_logs() -class LinChan_Adj(NekTestCase): - example_subdir = 'channel2D' - case_name = 'lin_chan_adj' +class LinCav_Adj(NekTestCase): + example_subdir = 'dfh_cav' + case_name = 'lin_dfh_cav_adj' def setUp(self): self.size_params = dict( ldim = '2', - lx1 = '10', - lxd = '15', + lx1 = '9', + lxd = '13', lx2 = 'lx1-2', - lelg = '50', + lelg = '500', lpelt = 'lelt', ) @@ -1006,11 +1006,11 @@ def setUp(self): def test_PnPn2_Serial(self): self.size_params['lx2'] = 'lx1-2' self.config_size() - self.build_nek(usr_file='lin_chan') + self.build_nek(usr_file='lin_dfh_cav') self.run_nek(step_limit=None) omega = self.get_value_from_log('Energy', column=-3, row=-1) - self.assertAlmostEqualDelayed(omega, target_val=-1.2337E-03, delta=2E-06, label='growth rate') + self.assertAlmostEqualDelayed(omega, target_val=-7.57304E-03, delta=1E-06, label='growth rate') self.assertDelayedFailures() @@ -1018,11 +1018,11 @@ def test_PnPn2_Serial(self): def test_PnPn2_Parallel(self): self.size_params['lx2'] = 'lx1-2' self.config_size() - self.build_nek(usr_file='lin_chan') + self.build_nek(usr_file='lin_dfh_cav') self.run_nek(step_limit=None) omega = self.get_value_from_log('Energy', column=-3, row=-1) - self.assertAlmostEqualDelayed(omega, target_val=-1.2337E-03, delta=2E-06, label='growth rate') + self.assertAlmostEqualDelayed(omega, target_val=-7.57304E-03, delta=1E-06, label='growth rate') self.assertDelayedFailures() diff --git a/short_tests/channel2D/lin_chan.usr b/short_tests/channel2D/lin_chan.usr deleted file mode 100644 index 9d184168a..000000000 --- a/short_tests/channel2D/lin_chan.usr +++ /dev/null @@ -1,206 +0,0 @@ -c----------------------------------------------------------------------- -c -c This is small 2D example of the linear direct and adjoint simulation -c of 2D Poiseuille flow, in which a fluid is moving laterally between -c two plates whose length and width is much greater than the distance -c separating them. Stability of this 2D parallel flow can be investigated -c analitcally by local analysis and we compare growth rate of the strongest -c mode with the stability calculation. -c -c----------------------------------------------------------------------- -c -c .Blasius I.C. and B.C. -c -c----------------------------------------------------------------------- - subroutine uservp (ix,iy,iz,ieg) - include 'SIZE' - include 'NEKUSE' ! UDIFF, UTRANS - - UDIFF =0. - UTRANS=0. - - return - end -c----------------------------------------------------------------------- - subroutine userf (ix,iy,iz,ieg) - include 'SIZE' - include 'NEKUSE' ! FFX, FFY, FFZ - - integer ix,iy,iz,ieg - - FFX = 0.0 - FFY = 0.0 - FFZ = 0.0 - - return - end -c----------------------------------------------------------------------- - subroutine userq (ix,iy,iz,ieg) - include 'SIZE' - include 'NEKUSE' ! QVOL - - QVOL = 0.0 - - return - end -c----------------------------------------------------------------------- - subroutine userchk - - include 'SIZE' ! NX1, NY1, NZ1, NELV, NIO - include 'INPUT' ! UPARAM - include 'TSTEP' ! ISTEP, IOSTEP, TIME, LASTEP - include 'SOLN' ! V[XYZ], V[XYZ]P, PRP, JP, VMULT - include 'MASS' ! BM1 - include 'ADJOINT' ! IFADJ - - integer n, nit - real Ek(2),timel(2), omega(2), domega - save EK, timel, omega - real vtmp(lx1*ly1*lz1*lelv,ldim) - character*132 restartf - -c set direct/adjoint mode and load restart field - if (ISTEP.eq.0) then - if (int(UPARAM(1)).eq.1) then - if (NIO.eq.0) write(*,*) 'Simulation in adjoint mode' - IFADJ = .TRUE. -c keep current base flow - call opcopy(vtmp(1,1),vtmp(1,2),vtmp(1,ndim),VX,VY,VZ) -c read the field - restartf = 'prtlin_chan_adj0.restart' - call load_fld(restartf) -c copy fileds - call opcopy(VXP,VYP,VZP,VX,VY,VZ) - n = NX2*NY2*NZ2*NELV - call copy(PRP,PR,n) -c put back base flow - call opcopy(VX,VY,VZ,vtmp(1,1),vtmp(1,2),vtmp(1,ndim)) - call rzero(PR,n) - else - if (NIO.eq.0) write(*,*) 'Simulation in direct mode' - IFADJ = .FALSE. -c keep current base flow - call opcopy(vtmp(1,1),vtmp(1,2),vtmp(1,ndim),VX,VY,VZ) -c read the field - restartf = 'prtlin_chan_dir0.restart' - call load_fld(restartf) -c copy fileds - call opcopy(VXP,VYP,VZP,VX,VY,VZ) - n = NX2*NY2*NZ2*NELV - call copy(PRP,PR,n) -c put back base flow - call opcopy(VX,VY,VZ,vtmp(1,1),vtmp(1,2),vtmp(1,ndim)) - call rzero(PR,n) - endif - endif - -c get energy - domega = 1.0 - nit = 10 - if (mod(ISTEP,nit).eq.0) then - n = NX1*NY1*NZ1*NELV - Ek(2) = Ek(1) - timel(2) = timel(1) - omega(2) = omega(1) - Ek(1) = 0.5*(glsc3(VXP,VXP,BM1,n)+glsc3(VYP,VYP,BM1,n)) - timel(1) = TIME -c get growthrate - if (Ek(2).gt.0.0.and.timel(1).gt.timel(2)) then - omega(1) = 0.5*log(Ek(1)/Ek(2))/(timel(1)-timel(2)) - domega = abs((omega(1)-omega(2))/omega(1)) - endif -c set logs - if (NIO.eq.0.and.ISTEP.gt.2*nit) - $ write(*,*) 'Energy ',Ek(1),omega(1),domega,TIME - endif - -c converged field - if (ISTEP.eq.NSTEPS) then -c write perturbation field - call outpost2(VXP,VYP,VZP,PRP,TP,0,'prt') - endif - - return - end -c----------------------------------------------------------------------- - subroutine userbc (ix,iy,iz,iside,eg) - - include 'SIZE' - include 'NEKUSE' ! UX, UY, UZ, TEMP, X, Y - -c velocity - UX = 0.0 - UY = 0.0 - UZ = 0.0 - - return - end -c----------------------------------------------------------------------- - subroutine useric (ix,iy,iz,ieg) - - include 'SIZE' - include 'NEKUSE' ! UX, UY, UZ, TEMP, Z - include 'SOLN' ! JP - - real amp, ran - -c velocity -c base flow - if (JP.eq.0) then - UX = (1.0-Y**2) - UY = 0.0 - UZ = 0.0 - else -c perturbation -c random distribution - amp = 1.0 - ran = 3.e4*(ieg+X*sin(Y)) - 1.5e3*ix*iy + .5e5*ix - ran = 1.e3*sin(ran) - ran = 1.e3*sin(ran) - ran = cos(ran) - UX = ran*amp - - ran = 2.3e4*(ieg+X*sin(Y)) + 2.3e3*ix*iy - 2.e5*ix - ran = 1.e3*sin(ran) - ran = 1.e3*sin(ran) - ran = cos(ran) - UY = ran*amp - - UZ = 0.0 - endif - - return - end -c----------------------------------------------------------------------- -c This routine to modify element vertices - subroutine usrdat - - include 'SIZE' - include 'TOTAL' - - integer n - real fact - - n = 8*nelv - - fact = 4.*atan(1.) - call cmult(xc,fact,n) - - x_min=glmin(xc,n) - y_min=glmin(yc,n) - x_max=glmax(xc,n) - y_max=glmax(yc,n) - - return - end -c----------------------------------------------------------------------- - subroutine usrdat2 - include 'SIZE' - - return - end -c----------------------------------------------------------------------- - subroutine usrdat3 - return - end -c----------------------------------------------------------------------- diff --git a/short_tests/channel2D/lin_chan_adj.re2 b/short_tests/channel2D/lin_chan_adj.re2 deleted file mode 100644 index 9988bedbb..000000000 Binary files a/short_tests/channel2D/lin_chan_adj.re2 and /dev/null differ diff --git a/short_tests/channel2D/lin_chan_dir.re2 b/short_tests/channel2D/lin_chan_dir.re2 deleted file mode 100644 index 9988bedbb..000000000 Binary files a/short_tests/channel2D/lin_chan_dir.re2 and /dev/null differ diff --git a/short_tests/channel2D/prtlin_chan_adj0.restart b/short_tests/channel2D/prtlin_chan_adj0.restart deleted file mode 100644 index 5faedce58..000000000 Binary files a/short_tests/channel2D/prtlin_chan_adj0.restart and /dev/null differ diff --git a/short_tests/channel2D/prtlin_chan_dir0.restart b/short_tests/channel2D/prtlin_chan_dir0.restart deleted file mode 100644 index 30bf62fd8..000000000 Binary files a/short_tests/channel2D/prtlin_chan_dir0.restart and /dev/null differ diff --git a/short_tests/channel2D/SIZE b/short_tests/dfh_cav/SIZE similarity index 90% rename from short_tests/channel2D/SIZE rename to short_tests/dfh_cav/SIZE index e71ee473c..898957a60 100644 --- a/short_tests/channel2D/SIZE +++ b/short_tests/dfh_cav/SIZE @@ -10,11 +10,11 @@ c ! BASIC parameter (ldim=2) ! domain dimension (2 or 3) - parameter (lx1=10) ! p-order (avoid uneven and values <6) - parameter (lxd=15) ! p-order for over-integration (dealiasing) + parameter (lx1=9) ! p-order (avoid uneven and values <6) + parameter (lxd=13) ! p-order for over-integration (dealiasing) parameter (lx2=lx1-2) ! p-order for pressure (lx1 or lx1-2) - parameter (lelg=50) ! max total number of elements + parameter (lelg=500) ! max total number of elements parameter (lpmin=1) ! min MPI ranks parameter (lpmax=1024) ! max MPI ranks parameter (ldimt=1) ! max auxiliary fields (temperature + scalars) diff --git a/short_tests/dfh_cav/baseflow.restart b/short_tests/dfh_cav/baseflow.restart new file mode 100644 index 000000000..425ef8095 Binary files /dev/null and b/short_tests/dfh_cav/baseflow.restart differ diff --git a/short_tests/dfh_cav/egvcavity_adj.restart b/short_tests/dfh_cav/egvcavity_adj.restart new file mode 100644 index 000000000..77a185b88 Binary files /dev/null and b/short_tests/dfh_cav/egvcavity_adj.restart differ diff --git a/short_tests/dfh_cav/egvcavity_dir.restart b/short_tests/dfh_cav/egvcavity_dir.restart new file mode 100644 index 000000000..3331501a6 Binary files /dev/null and b/short_tests/dfh_cav/egvcavity_dir.restart differ diff --git a/short_tests/dfh_cav/lin_dfh_cav.usr b/short_tests/dfh_cav/lin_dfh_cav.usr new file mode 100644 index 000000000..545010eed --- /dev/null +++ b/short_tests/dfh_cav/lin_dfh_cav.usr @@ -0,0 +1,379 @@ +c----------------------------------------------------------------------- +c +c This is small 2D example of the linear direct and adjoint simulation +c of differentially heated square cavity. The linear stability of natural +c convection is tested and the grow rate of one of the modes is compared +c with the result of the linear stability analysis. +c This simulation tests temperature adjoint. +c +c----------------------------------------------------------------------- +c +c .Blasius I.C. and B.C. +c +c----------------------------------------------------------------------- + subroutine uservp (ix,iy,iz,ieg) + include 'SIZE' + include 'NEKUSE' ! UDIFF, UTRANS + include 'INPUT' ! UPARAM + include 'ADJOINT' ! IFADJ + include 'TSTEP' ! TSTEP + +c argument list + integer ix, iy, iz, ieg + +c local variables + real Pr_,Ra_ + +c decide the non dimensionalization of the equations + Pr_ = abs(UPARAM(2)) + Ra_ = abs(UPARAM(3)) +c set the coefficients +c adjoint problem + if (IFADJ) then + if (IFIELD.eq.1) then ! momentum equations + UTRANS = 1.0 + UDIFF = Pr_/sqrt(Ra_) + elseif (IFIELD.eq.2) then ! temperature equation + UTRANS = 1.0 + UDIFF = 1./sqrt(Ra_) + endif +c direct problem + else + if (IFIELD.eq.1) then ! momentum equations + UTRANS = 1./Pr_ + UDIFF = 1./sqrt(Ra_) + elseif (IFIELD.eq.2) then ! temperature equation + UTRANS = 1.0 + UDIFF = 1./sqrt(Ra_) + endif + endif + + return + end +c----------------------------------------------------------------------- + subroutine userf (ix,iy,iz,ieg) + include 'SIZE' + include 'NEKUSE' ! FF[XYZ] + include 'PARALLEL' ! GLLEL + include 'SOLN' ! VTRANS,TP + include 'INPUT' ! IF3D + include 'ADJOINT' ! IFADJ, G_ADJ + +c argument list + integer ix,iy,iz,ieg + +c local variable + integer ip,iel + real rtmp + +c local element number + iel=GLLEL(ieg) +c initialize the forcing to zero + FFX = 0.0 + FFY = 0.0 + if (IF3D) FFZ = 0.0 +c compute the corresponding index in the pertubation arrays + ip=ix+NX1*(iy-1+NY1*(iz-1+NZ1*(iel-1))) +c forcing for the direct run is the buoyancy expressed in terms of +c Boussinesq approximation + if (.not.IFADJ) then + rtmp = TP(ip,1,1)/VTRANS(ix,iy,iz,iel,1) + FFX = FFX + G_ADJ(1)*rtmp + FFY = FFY + G_ADJ(2)*rtmp + if (IF3D) FFZ = FFZ + G_ADJ(3)*rtmp + endif + + return + end +c----------------------------------------------------------------------- + subroutine userq (ix,iy,iz,ieg) + include 'SIZE' + include 'NEKUSE' ! QVOL + include 'PARALLEL' ! GLLEL + include 'SOLN' ! VTRANS,V[XYZ]P + include 'INPUT' ! IF3D + include 'ADJOINT' ! IFADJ, G_ADJ, BETA_B + +c argument list + integer ix,iy,iz,ieg + +c local variable + integer ip,iel + real rtmp + +c local element number + iel = GLLEL(ieg) +c initialize the forcing to zero + QVOL = 0.0 +c compute the corresponding index in the pertubation arrays + ip = ix+NX1*(iy-1+NY1*(iz-1+NZ1*(iel-1))) +c forcing for the adjoint temperature equation that depends on the +c approximation used for the buoyancy (in theis case Boussinesq) + if (IFADJ) then + QVOL = (VXP(ip,1)*G_ADJ(1)+VYP(ip,1)*G_ADJ(2) + $ +VZP(ip,1)*G_ADJ(3))*BETA_B(ip) !const g/|g|.u + endif + + return + end +c----------------------------------------------------------------------- + subroutine userchk + + include 'SIZE' ! NX1, NY1, NZ1, NELV, NIO + include 'INPUT' ! UPARAM, CPFLD + include 'TSTEP' ! ISTEP, IOSTEP, TIME, LASTEP + include 'SOLN' ! V[XYZ], T, V[XYZ]P, TP, PRP + include 'MASS' ! BM1 + include 'ADJOINT' ! IFADJ, G_ADJ, BETA_B, DTD[XYZ] + +c local variables + integer n, nit + real Pr_,Ra_ + real Ek(2),timel(2), omega(2), domega + save EK, timel, omega + real vtmp(lx1*ly1*lz1*lelt,ldim),ttmp(lx1*ly1*lz1*lelt) + character*132 restartf + real cht_glsc2_wt + + if (ISTEP.eq.0) then +c set direct/adjoint mode and load restart field +c keep current base flow + call opcopy(vtmp(1,1),vtmp(1,2),vtmp(1,ndim),VX,VY,VZ) + n = NX1*NY1*NZ1*NELV + call copy(ttmp,T,n) + + if (int(UPARAM(1)).eq.1) then + if (NIO.eq.0) write(*,*) 'Simulation in adjoint mode' + IFADJ = .TRUE. +c restart file name + restartf = 'egvcavity_adj.restart' + else + if (NIO.eq.0) write(*,*) 'Simulation in direct mode' + IFADJ = .FALSE. +c restart file name + restartf = 'egvcavity_dir.restart' + endif + +c read the field + call load_fld(restartf) +c copy fileds + call opcopy(VXP,VYP,VZP,VX,VY,VZ) + n = NX1*NY1*NZ1*NELV + call copy(TP,t,n) + n = NX2*NY2*NZ2*NELV + call copy(PRP,PR,n) +c put back base flow + call opcopy(VX,VY,VZ,vtmp(1,1),vtmp(1,2),vtmp(1,ndim)) + n = NX1*NY1*NZ1*NELV + call copy(T,ttmp,n) + call rzero(PR,n) + +c decide the non dimensionalization of the equations + Pr_ = abs(UPARAM(2)) + Ra_ = abs(UPARAM(3)) + + +c temperature base flow gradient, needed for the extra advection +c terms in the adjoint momentum equations + if (IFADJ) then + call gradm1(DTDX,DTDY,DTDZ,T) +c coefficient in front of the buoyancy term, depends on the +c nondimensionalization + call rone(BETA_B,n) + call cmult(BETA_B,Pr_,n) + endif + +c gravity + G_ADJ(1) = 0.0 !x direction + G_ADJ(2) = 1.0 !y direction + G_ADJ(3) = 0.0 !z direction + +c set fluid properties + if (IFADJ) then + CPFLD(1,1)=Pr_/sqrt(Ra_) + CPFLD(1,2)=1.0 + + CPFLD(2,1)=1.0/sqrt(Ra_) + CPFLD(2,2)=1.0 + else + CPFLD(1,1)=1.0/sqrt(Ra_) + CPFLD(1,2)=1.0/Pr_ + + CPFLD(2,1)=1.0/sqrt(Ra_) + CPFLD(2,2)=1.0 + endif + endif + +c get energy + domega = 1.0 + nit = 10!00 + if (mod(ISTEP,nit).eq.0) then + n = NX1*NY1*NZ1*NELV + Ek(2) = Ek(1) + timel(2) = timel(1) + omega(2) = omega(1) +! Ek(1) = 0.5*(glsc3(VXP,VXP,BM1,n)+glsc3(VYP,VYP,BM1,n)) + Ek(1) = cht_glsc2_wt(VXP,VYP,VZP,TP,VXP,VYP,VZP,TP,BM1) + timel(1) = TIME +c get growth rate + if (Ek(2).gt.0.0.and.timel(1).gt.timel(2)) then + omega(1) = 0.5*log(Ek(1)/Ek(2))/(timel(1)-timel(2)) + domega = abs((omega(1)-omega(2))/omega(1)) + endif +c set logs + if (NIO.eq.0.and.ISTEP.gt.2*nit) + $ write(*,*) 'Energy ',Ek(1),omega(1),domega,TIME + endif + +c converged field + if (ISTEP.eq.NSTEPS) then +c write perturbation field + call outpost2(VXP,VYP,VZP,PRP,TP,0,'prt') + endif + + return + end +c----------------------------------------------------------------------- + subroutine userbc (ix,iy,iz,iside,eg) + + include 'SIZE' + include 'NEKUSE' ! UX, UY, UZ, TEMP, X, Y + +c velocity + UX = 0.0 + UY = 0.0 + UZ = 0.0 + +c temperature + TEMP = 0.0 + + return + end +c----------------------------------------------------------------------- + subroutine useric (ix,iy,iz,ieg) + + include 'SIZE' + include 'NEKUSE' ! UX, UY, UZ, TEMP, Z + + real amp, ran + +c velocity random distribution + amp = 0.001 + + ran = 3.e4*(ieg+X*sin(Y)) - 1.5e3*ix*iy + .5e5*ix + ran = 1.e3*sin(ran) + ran = 1.e3*sin(ran) + ran = cos(ran) + UX = ran*amp + + ran = 2.3e4*(ieg+X*sin(Y)) + 2.3e3*ix*iy - 2.e5*ix + ran = 1.e3*sin(ran) + ran = 1.e3*sin(ran) + ran = cos(ran) + UY = ran*amp + + UZ = 0.0 + +c temperature random distribution + ran = 2.e4*(ieg+x*sin(y)) + 1.e3*ix*iy + 1.e5*ix + ran = 1.e3*sin(ran) + ran = 1.e3*sin(ran) + ran = cos(ran) + TEMP= ran*amp + + return + end +c----------------------------------------------------------------------- +c This routine to modify element vertices + subroutine usrdat + include 'SIZE' + + return + end +c----------------------------------------------------------------------- + subroutine usrdat2 + include 'SIZE' + + return + end +c----------------------------------------------------------------------- + subroutine usrdat3 + return + end +c----------------------------------------------------------------------- +c global scalar (L2 norm) + real function cht_glsc2_wt (b1,b2,b3,b4,x1,x2,x3,x4,wt) + implicit none + + include 'SIZE' ! N[XYZ]1, NEL[VT] + include 'INPUT' ! IFFLOW, IFHEAT, IF3D + include 'MASS' ! VOLVM1, VOLTM1 + +c argument list + real b1(1),b2(1),b3(1),b4(1),x1(1),x2(1),x3(1),x4(1),wt(1) + +c local variables + integer ntotv, ntott + real sum, f1, f2 + integer i + real CHCST_SC,CHCFF_V,CHCFF_T + +c functions + real glsum +c----------------------------------------------------------------------- + ntotv = NX1*NY1*NZ1*NELV + ntott = NX1*NY1*NZ1*NELT + +c default values + CHCST_SC = 3.36558d0 + CHCFF_V = 0.5 + CHCFF_T = 0.5 + +c scaling factor velocity vs temperature + f1=CHCFF_V/VOLVM1 + f2=CHCFF_T*CHCST_SC/VOLTM1 + + sum = 0. + if (IFFLOW) then !if vel + if (IFHEAT) then !if temp & vel + if (IF3D) then + do i=1,ntotv + sum = sum + wt(i)*(f1*(b1(i)*x1(i)+b2(i)*x2(i) + & +b3(i)*x3(i))+f2*b4(i)*x4(i)) + end do + else + do i=1,ntotv + sum =sum + wt(i)*(f1*(b1(i)*x1(i)+b2(i)*x2(i)) + & +f2*b4(i)*x4(i)) + end do + end if + if (ntott.gt.ntotv) then !conjugate heat transfer + do i=ntotv+1,ntott + sum = sum + wt(i)*f2*b4(i)*x4(i) + end do + end if + else !just vel + if (IF3D) then + do i=1,ntotv + sum = sum + wt(i)*f1*(b1(i)*x1(i)+ + $ b2(i)*x2(i)+b3(i)*x3(i)) + end do + else + do i=1,ntotv + sum = sum + wt(i)*f1*(b1(i)*x1(i)+b2(i)*x2(i)) + end do + end if + end if + else !just temp + if (IFHEAT) then + do i=1,ntott + sum = sum + wt(i)*(f2*b4(i)*x4(i)) + end do + end if + end if + + cht_glsc2_wt = glsum(sum,1) + + return + end +c----------------------------------------------------------------------- diff --git a/short_tests/channel2D/lin_chan_dir.par b/short_tests/dfh_cav/lin_dfh_cav_adj.par similarity index 50% rename from short_tests/channel2D/lin_chan_dir.par rename to short_tests/dfh_cav/lin_dfh_cav_adj.par index b10b88ee3..89fa51526 100644 --- a/short_tests/channel2D/lin_chan_dir.par +++ b/short_tests/dfh_cav/lin_dfh_cav_adj.par @@ -1,9 +1,12 @@ [GENERAL] +startFrom = baseflow.restart stopat = numSteps #endTime numsteps = 1000 -perturbationmodes = -1 -userparam01 = 0 -dt = -1.0e-02 +perturbationmodes = 1 +userparam01 = 1 # 0-direct; 1-adjoint +userparam02 = 0.01500 # PRANDTL number +userparam03 = 45000.0 # RAYLEIGH number +dt = -5.0e-03 timestepper = bdf #char #steady torder = 2 writecontrol = timeStep #runTime @@ -14,7 +17,8 @@ filterweight = 1e-02 [PROBLEMTYPE] perturbations = yes -solvebaseflow = no +solveBaseflow = no +variableProperties = yes [PRESSURE] preconditioner = semg @@ -25,4 +29,10 @@ residualproj = no residualtol = 1e-08 residualproj = no density = 1 -viscosity = -2000.0 +viscosity = 1.0 + +[TEMPERATURE] # temperature with Hmholtz +conjugateHeatTransfer = no +rhoCp = 1.0 +conductivity = 0.5 +residualTol = 1e-08 diff --git a/short_tests/dfh_cav/lin_dfh_cav_adj.re2 b/short_tests/dfh_cav/lin_dfh_cav_adj.re2 new file mode 100644 index 000000000..7bf533964 Binary files /dev/null and b/short_tests/dfh_cav/lin_dfh_cav_adj.re2 differ diff --git a/short_tests/channel2D/lin_chan_adj.par b/short_tests/dfh_cav/lin_dfh_cav_dir.par similarity index 50% rename from short_tests/channel2D/lin_chan_adj.par rename to short_tests/dfh_cav/lin_dfh_cav_dir.par index 6d0cd461c..b8a7ef773 100644 --- a/short_tests/channel2D/lin_chan_adj.par +++ b/short_tests/dfh_cav/lin_dfh_cav_dir.par @@ -1,9 +1,12 @@ [GENERAL] +startFrom = baseflow.restart stopat = numSteps #endTime numsteps = 1000 -perturbationmodes = -1 -userparam01 = 1 -dt = -1.0e-02 +perturbationmodes = 1 +userparam01 = 0 # 0-direct; 1-adjoint +userparam02 = 0.01500 # PRANDTL number +userparam03 = 45000.0 # RAYLEIGH number +dt = -5.0e-03 timestepper = bdf #char #steady torder = 2 writecontrol = timeStep #runTime @@ -14,7 +17,8 @@ filterweight = 1e-02 [PROBLEMTYPE] perturbations = yes -solvebaseflow = no +solveBaseflow = no +variableProperties = yes [PRESSURE] preconditioner = semg @@ -25,4 +29,10 @@ residualproj = no residualtol = 1e-08 residualproj = no density = 1 -viscosity = -2000.0 +viscosity = 1.0 + +[TEMPERATURE] # temperature with Hmholtz +conjugateHeatTransfer = no +rhoCp = 1.0 +conductivity = 0.5 +residualTol = 1e-08 \ No newline at end of file diff --git a/short_tests/dfh_cav/lin_dfh_cav_dir.re2 b/short_tests/dfh_cav/lin_dfh_cav_dir.re2 new file mode 100644 index 000000000..7bf533964 Binary files /dev/null and b/short_tests/dfh_cav/lin_dfh_cav_dir.re2 differ