diff --git a/standalone/fuel_interp_test_main.F b/standalone/fuel_interp_test_main.F index f2856c51..39baf4df 100644 --- a/standalone/fuel_interp_test_main.F +++ b/standalone/fuel_interp_test_main.F @@ -2,24 +2,20 @@ program test use module_fr_sfire_core implicit none real :: time_now -integer :: ims,ime,jms,jme,icl,jcl,i,j,num_tests +integer :: ims,ime,jms,jme,icl,jcl,i,j,num_tests,flag real,dimension(3,3)::tff,lff,lfn,tign,real_tff,real_lff real,dimension(3,3)::tff0,lff0,lfn0,tign0,real_tff0,real_lff0 real,dimension(3,3)::tff1,lff1,lfn1,tign1,real_tff1,real_lff1 real,dimension(3,3)::tff2,lff2,lfn2,tign2,real_tff2,real_lff2 real,dimension(3,3)::tff3,lff3,lfn3,tign3,real_tff3,real_lff3 -real,dimension(3,3)::tff4,lff4,lfn4,tign4 -real,dimension(3,3)::tff5,lff5,lfn5,tign5 -real,dimension(3,3)::tff6,lff6,lfn6,tign6 -real,dimension(3,3)::tff7,lff7,lfn7,tign7 real,dimension(3,3)::tff8,lff8,lfn8,tign8 real,dimension(3,5)::tff9,lff9,lfn9,tign9 real, dimension(5,5)::d real,dimension(2,2)::result -real::test_err,error,test1_err,test2_err,test3_err,test4_err,test5_err +real::test_err,error real::lfn00,lfn01,lfn10,lfn11,tign00,tign01,tign10,tign11,fuel_time_cell -real:: res,Result1,Result2,Result3,Result4,Result5,result6 -real:: res1,res2,res3,res4,res5 +real:: Result1,Result2,Result3,Result4,Result5,result6 +real:: res,res1,res2,res3,res4! time_now=3. 101 format(a/(3f7.3)) ims=1 @@ -368,19 +364,33 @@ program test end do end do +if (res1-result(1,1)>error.or.res2-result(1,2)>error.or.res3-result(2,1)>error.or.res4-result(2,2)>error) then + print*,'Test #10 is not working, one of errors ' + print*,'For cell (1,1)',res1-result(1,1) + print*,'For cell (1,2)',res2-result(1,2) + print*,'For cell (2,1)',res3-result(2,1) + print*,'For cell (2,2)',res4-result(2,2) + print*,'is greater than 0.0001',error +else +num_tests=num_tests+1; +endif print*,"" print*,"" !!!! Test #2 -write(*,*)'MAIN TEST #2' +write(*,*)'MAIN TEST #2 (#11)' print*,"Check for consistency, for 3 inner cells with the same results" time_now=2. ims=1 ime=5 jms=1 jme=5 - +res1=2.9752314E-02 +res2=1.000 +res3=2.9752314E-02 +res4=1.000 +flag=1 data lfn9/-1., -1.,-1., & -1., -1.,0., & 0., 0.,0., & @@ -396,19 +406,7 @@ program test do icl=2,4 write(*,*)"icl,jcl",icl,jcl call tign_lfn_interpolation(time_now,icl,jcl,ims,ime,jms,jme, & - tign9,lfn9,tff,lff) - -d=1e20 -d(1:5:2,1:5:2)=tign9(icl-1:icl+1,1:3) -d(2:4,2:4)=tff -!write(*,*)'tign and tff' -!write(*,'(5f7.3)')d -!write(*,*)' ' -d(1:5:2,1:5:2)=lfn9(icl-1:icl+1,1:3) -d(2:4,2:4)=lff -!write(*,*)'lfn and lff' -!write(*,'(5f7.3)')d -!write(*,*)' ' + tign9,lfn9,tff9,lff9) fuel_time_cell= 8.235294 ; write(*,*)"Calculation of fuel_frac over 4 subcells" @@ -417,15 +415,29 @@ program test do j=1,2 -result6=fuel_left_cell_3( & - lff(i,j),lff(i,j+1),lff(i+1,j),lff(i+1,j+1), & - tff(i,j),tff(i,j+1),tff(i+1,j),tff(i+1,j+1),& +result(i,j)=fuel_left_cell_3( & + lff9(i,j),lff9(i,j+1),lff9(i+1,j),lff9(i+1,j+1), & + tff9(i,j),tff9(i,j+1),tff9(i+1,j),tff9(i+1,j+1),& time_now, fuel_time_cell) -write(*,*)'result after step i and j',result,i,j +write(*,*)'result after step i and j',result(i,j),i,j end do end do + +if (res1-result(1,1)>error.or.res2-result(1,2)>error.or.res3-result(2,1)>error.or.res4-result(2,2)>error) then + print*,'Test #10 is not working, one of errors ' + print*,'For cell (1,1)',res1-result(1,1) + print*,'For cell (1,2)',res2-result(1,2) + print*,'For cell (2,1)',res3-result(2,1) + print*,'For cell (2,2)',res4-result(2,2) + print*,'is greater than 0.0001',error + flag=0 +end if + end do +if (flag.eq.1) then + num_tests=num_tests+1; +endif print*,"If the results for each correspondent subcell are equal in all cells" print*,"then the code is consistent" diff --git a/wrfv2_fire/phys/module_fr_sfire_core.F b/wrfv2_fire/phys/module_fr_sfire_core.F index 4ce3f497..52d63e6c 100644 --- a/wrfv2_fire/phys/module_fr_sfire_core.F +++ b/wrfv2_fire/phys/module_fr_sfire_core.F @@ -263,7 +263,6 @@ subroutine fuel_left( & real,dimension(3,3)::tff,lff ! help variables instead of arrays fuel_left_f and fire_area_f real::lffij,lffi1j,lffij1,lffi1j1,tifij,tifi1j,tifij1,tifi1j1,tx,ty,txx,tyy -! variables for calculation instead of lff(i,j) and tif(i,j)is lffij,tifij etc..#define IFCELLS (ite-its+1)*fuel_left_irl #define JFCELLS (jte-jts+1)*fuel_left_jrl character(len=128)::msg integer::m,omp_get_thread_num @@ -306,25 +305,14 @@ subroutine fuel_left( & ite1=min(ite,ifde-1) jts1=max(jts,jfds+1) jte1=min(jte,jfde-1) -!write(*,*)"fuel_left_method",fuel_left_method do icl=its1,ite1 do jcl=jts1,jte1 helpsum1=0 helpsum2=0 - - ! Loop over subcells in cell #(icl,jcl) -! write(*,*)"ifds,ifde",ifds,ifde -! write(*,*)"jfds,jfde",jfds,jfde -! write(*,*)"its,ite,jts,jte",its,ite,jts,jte -! write(*,*)"icl,jcl",icl,jcl - - call tign_lfn_interpolation(time_now,icl,jcl,ims,ime,jms,jme, & tign,lfn,tff,lff) -!endif - do isubcl=1,ir do jsubcl=1,jr !we use 4 points of each subcell in fuel_left_cell_1 @@ -340,7 +328,6 @@ subroutine fuel_left( & lff(isubcl,jsubcl),lff(isubcl,jsubcl+1),lff(isubcl+1,jsubcl),lff(isubcl+1,jsubcl+1), & tff(isubcl,jsubcl),tff(isubcl,jsubcl+1),tff(isubcl+1,jsubcl),tff(isubcl+1,jsubcl+1), & time_now, fuel_time(icl,jcl)) -! dont forget to change fire_area_ff here else call crash('fuel_left: unknown fuel_left_method') endif @@ -504,24 +491,26 @@ end subroutine tign_lfn_interpolation subroutine tign_line_interp(tign1,tign2,lfn1,lfn2,lfn_subcl,tign_subcl,time_now) +! Description: subroutine that is calculating tign of subcell that lies on the +! line between 2 cells, lfn_subcl is interpolated linearly look at tign_lfn_interpolation real,intent(in) :: tign1,tign2,lfn1,lfn2,lfn_subcl,time_now real,intent(out) :: tign_subcl real :: a,b,c,err err=0.1 -! write(*,*)"lfn1,lfn2,tign1,tign2,lfn_subcl,timenow",lfn1,lfn2,tign1,tign2,lfn_subcl,time_now if((lfn1.le.0).AND.(lfn2.le.0)) then +! both lfn1,2<0, so we do linear interpolation tign_subcl=0.5*(tign1+tign2) elseif((lfn1.gt.0).AND.(lfn2.gt.0))then +! both lfn1,2>0, so we do linear interpolation, so it is equal to time_now if ((abs(tign1-time_now).gt.err).or.(abs(tign2-time_now).gt.err)) then - write(*,*)"lfn1,lfn2,tign1,tign2,timenow",lfn1,lfn2,tign1,tign2,time_now call crash('line_interp: when lfn1(2)>0 their tign1(2) should = time_now') else tign_subcl=time_now; endif -elseif(lfn_subcl.gt.0) then +elseif(lfn_subcl.gt.0) then !So the subcell is not on fire if ((abs(tign1-time_now).gt.err).OR.(abs(tign2-time_now).gt.err)) then call crash('tign_line_interp one of tign1,2 should be equal time_now') else @@ -555,7 +544,9 @@ end subroutine tign_line_interp subroutine tign_four_pnts_interp(tign1,tign2,tign3,tign4,lfn1,lfn2, & lfn3,lfn4,lfn_subcl,tign_subcl,time_now) - +! Description: subroutine that is calculating tign of subcell that lies in the +! middle of the square with 4 cells in the endpoints, lfn_subcl is interpolated linearly +! look at tign_lfn_interpolation real,intent(in) :: tign1,tign2,tign3,tign4 real,intent(in) :: lfn1,lfn2,lfn3,lfn4,lfn_subcl,time_now real,intent(out) :: tign_subcl @@ -564,14 +555,16 @@ subroutine tign_four_pnts_interp(tign1,tign2,tign3,tign4,lfn1,lfn2, & err=0.0001 if((lfn1.le.0).AND.(lfn2.le.0).AND.(lfn3.le.0).AND.(lfn4.le.0)) then +! all 4 cells:i:1,4 lfn(i)<0, so we do linear interpolation tign_subcl=0.25*(tign1+tign2+tign3+tign4) elseif((lfn1.gt.0).AND.(lfn2.gt.0).AND.(lfn3.gt.0).AND.(lfn4.gt.0))then +! all 4 cells:i:1,4 lfn(i)>0, not on fire so tign_subcl=time_now !if ((abs(tign1-time_now).gt.err).OR.(abs(tign2-time_now).gt.err).OR.(abs(tign3-time_now).gt.err).OR.(abs(tign4-time_now).gt.err)) then ! call crash('tign_four_pnts_interp: when lfn1(2,3,4)>0 their tign1(2,3,4) should = time_now') !else tign_subcl=time_now; !endif -elseif(lfn_subcl.gt.0) then +elseif(lfn_subcl.gt.0) then ! so the subcell is not on fire ! if ((abs(tign1-time_now).gt.err).OR.(abs(tign2-time_now).gt.err).OR.(abs(tign3-time_now).gt.err).OR.(abs(tign4-time_now).gt.err)) then ! call crash('tign_line_interp one of tign1(2,3,4) should be equal time_now') ! else @@ -583,6 +576,7 @@ subroutine tign_four_pnts_interp(tign1,tign2,tign3,tign4,lfn1,lfn2, & ! tign_subcl~=c*lfn_subcl+time_now; a=0; b=0; +!Calculating weights if (lfn1.le.0) then ! if (tign1.gt.time_now) ! call crash('tign_four_pnts_interp tign1 should be less then time_now'); @@ -936,14 +930,9 @@ real function fuel_left_cell_3( & external dggglm !executable statements -! a very crude approximation - replace by a better code -! call check_mesh_2dim(ids,ide+1,jds,jde+1,ims,ime,jms,jme) dx=1 dy=1 -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!comment - changed tign-time_now to time_now-tign -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! t00=time_now-tign00 if(lfn00>=0. .or. t00<0.)t00=0. t01=time_now-tign01 @@ -954,14 +943,10 @@ real function fuel_left_cell_3( & if(lfn11>=0. .or. t11<0.)t11=0. !*** case0 Do nothing if ( lfn00>=0 .and. lfn10>=0 .and. lfn01>=0 .and. lfn11>=0 ) then -!print*,"Case 0" out = 1.0 ! fuel_left, no burning -!*** case4 all four coners are burning +!*** case4 all four corners are burning else if (lfn00<=0 .and. lfn10<=0 .and. lfn01<=0 .and. lfn11<=0) then -!!!!!!!!!!! - -!print*,"Case 4" -! All burning +! All burning, use Finite Differences ! T=u(1)*x+u(2)*y+u(3) ! t(0,0)=tign(1,1) ! t(0,fd(2))=tign(1,2) @@ -977,41 +962,28 @@ real function fuel_left_cell_3( & ! u(3)=t(h/2,h/2)-h/2(dt/dx+dt/dy) tign_middle=(t00+t01+t10+t11)/4 -! write(*,*)"tign_middle",tign_middle -! write(*,*)"t00,t01,t10,t11",t00,t01,t10,t11 ! since mesh_size is 1 we replace fd(1) and fd(2) by 1 dt_dx=(t10-t00+t11-t01)/2 dt_dy=(t01-t00+t11-t10)/2 -! write(*,*)"dt_dx,dt_dy",dt_dx,dt_dy -! write(*,*)"dx,dy",dx,dy - u(1)=dt_dx u(2)=dt_dy u(3)=tign_middle-(dt_dx+dt_dy)/2 -!Write(*,*)"coefficients of the plane u=",u(1),u(2),u(3) - - ! integrate - u(1)=-u(1)/fuel_time_cell - u(2)=-u(2)/fuel_time_cell - u(3)=-u(3)/fuel_time_cell -!write(*,*)"u/fuel_time_cell",u(1),u(2),u(3) - -!write(*,*)"intexp(u(1)),intexp(u(2))",intexp(u(1)*dx),intexp(u(2)*dy) - !fuel_burn(i,j)=1-exp(u(3))*intexp(u(1)*dx)*intexp(u(2)*dy) + ! integrate + u(1)=-u(1)/fuel_time_cell + u(2)=-u(2)/fuel_time_cell + u(3)=-u(3)/fuel_time_cell s1=u(1) s2=u(2) out=1-exp(u(3))*intexp(s1)*intexp(s2) -! print *,intexp(s1),intexp(s2),out if ( out<0 .or. out>1.0 ) then print *,'case4, out should be between 0 and 1' end if -!*** case 1,2,3 +!*** case 1,2,3- other cases else -!print*,"Case 123" ! set xx, yy for the coner points ! move these values out of i and j loop to speed up xx(1) = -0.5 @@ -1065,24 +1037,19 @@ real function fuel_left_cell_3( & xylist(npoint+1,2)=xylist(1,2) - - !print *,'after LS in case3'pproximation of the plane for lfn using finite - !differences +!approximation of the plane for lfn using finite differences, look at "case 4" explanations ! approximate L(x,y)=u(1)*x+u(2)*y+u(3) lfn_middle=(lfn00+lfn01+lfn10+lfn11)/4 -! print *,"lfn_middle",lfn_middle dt_dx=(lfn10-lfn00+lfn11-lfn01)/2 dt_dy=(lfn01-lfn00+lfn11-lfn10)/2 -!print *,"dt_dx,dt_dy",dt_dx,dt_dy u(1)=dt_dx u(2)=dt_dy u(3)=lfn_middle-(dt_dx+dt_dy)/2 ! finding the coefficient c, reminder we work over one subcell only ! T(x,y)=c*L(x,y)+time_now -!write(*,*)"vector u before c",u(1),u(2),u(3) a=0 b=0 - +!Calculating weights if (lfn00 <= 0) then a=a+lfn00*lfn00 if (t00 < 0) then @@ -1090,8 +1057,6 @@ real function fuel_left_cell_3( & else b=b+t00*lfn00 end if - ! write(*,*)"t00,time_now,lfn00",t00,time_now,lfn00 - ! write(*,*)"a,b",a,b end if @@ -1102,7 +1067,6 @@ real function fuel_left_cell_3( & else b=b+t01*lfn01 end if - ! write(*,*)"a,b",a,b end if @@ -1113,7 +1077,6 @@ real function fuel_left_cell_3( & else b=b+t10*lfn10 end if - ! write(*,*)"a,b",a,b end if if (lfn11<=0) then @@ -1123,7 +1086,6 @@ real function fuel_left_cell_3( & else b=b+t11*lfn11 end if - ! write(*,*)"a,b",a,b end if @@ -1136,12 +1098,6 @@ real function fuel_left_cell_3( & u(2)=u(2)*c u(3)=u(3)*c -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - !print *,'vecX from LS',vecX - !print *,'tign inputed',tign00,tign10,tign11,tign01 - -!write(*,*)"vector u",u(1),u(2),u(3) ! rotate to gradient on x only nr = sqrt(u(1)**2+u(2)**2) if(.not.nr.gt.eps)then