demona: D:/demona/2.0/manual/doc/2-31set/problems/ldc

00001 module model
00002 !module modelSet  ! use to differentiate between set and no set models.
00003 
00004 
00005 use global
00006 
00007 use myFunctions
00008 
00009 implicit none
00010 
00011 
00012 contains 
00013 
00014 
00015 function discrete(a)
00016 
00017 
00018 real(8) discrete(globalCurrentSetLength)
00019 
00020 integer i,j,t
00021 
00022 integer a
00023 
00024 integer myI,myJ,myT,myBC,myUnk
00025 
00026 integer myP,myD
00027 
00028 integer,pointer :: IM(:,:,:,:)
00029 
00030 integer,pointer :: IM2(:,:,:,:)
00031 
00032 integer,pointer :: maxI,maxJ,maxT
00033 
00034 integer,pointer :: iVec(:),jVec(:),bcVec(:),tVec(:)
00035 
00036 integer,pointer :: endIndex(:)
00037 
00038 integer modelPhysics
00039 integer modelDomain
00040 
00041 real(8), allocatable :: RHS(:)
00042 
00043 integer,pointer :: dof
00044 
00045 integer,pointer :: nx,ny
00046 
00047 !integer,pointer :: dof
00048 
00049 
00050 ! model specific definitions
00051 
00052 real(8) PC,PW,PE,PS,PN
00053 real(8) UC,UW,UE,US,UN
00054 real(8) VC,VW,VE,VS,VN
00055 
00056 real(8) PEE,PWW,PNN,PSS
00057 
00058 real(8) UEE,UWW
00059 real(8) UEEE,UWWW
00060 real(8) UNN,USS
00061 real(8) UNNN,USSS
00062 real(8) VEE,VWW
00063 real(8) VEEE,VWWW
00064 real(8) VNN,VSS
00065 real(8) VNNN,VSSS
00066 
00067 real(8) rhoC,rhoW,rhoE,rhoS,rhoN
00068 real(8) muC,muW,muE,muS,muN
00069 
00070 real(8) rho,mu
00071 
00072 
00073 real(8) PCPTO
00074 
00075 real(8) UCPTO,UCRTO
00076 
00077 real(8) VCPTO,VCRTO
00078 
00079 
00080 real(8) dPdX,dPdY
00081 real(8) dUdX,dUdY
00082 real(8) dVdX,dVdY
00083 real(8) d2UdX2,d2UdY2
00084 real(8) d2VdX2,d2VdY2
00085 
00086 real(8) dPdTau
00087 
00088 real(8) dUdTau,dUdT
00089 
00090 real(8) dVdTau,dVdT
00091 
00092 
00093 
00094 real(8) UdUdX,VdUdY
00095 real(8) UdVdX,VdVdY
00096 
00097 
00098 ! finite volume discretizations
00099 real(8) dUdXe,dUdXw,dUdYn,dUdYs
00100 real(8) dVdXe,dVdXw,dVdYn,dVdYs
00101 
00102 real(8) U_fd_e,U_fd_w,U_fd_n,U_fd_s
00103 real(8) V_fd_e,V_fd_w,V_fd_n,V_fd_s
00104 
00105 real(8) U_fd,V_fd
00106 
00107 real(8) U_fc_e,U_fc_w,U_fc_n,U_fc_s
00108 real(8) V_fc_e,V_fc_w,V_fc_n,V_fc_s
00109 
00110 real(8) U_fc,V_fc
00111 
00112 real(8) Ufw,Ufe,Ufn,Ufs
00113 real(8) Vfw,Vfe,Vfn,Vfs
00114 real(8) rhofw,rhofe,rhofn,rhofs
00115 
00116 
00117 real(8) Px,Py
00118 
00119 real(8) upw1,upw2,down,value
00120 
00121 
00122 ! MLN parameters
00123 
00124 real(8) au,Ap
00125 
00126 
00127 
00128 !real(8) dtau
00129 
00130 real(8) upu,upv
00131 
00132 real(8) myU,myV
00133 
00134 real(8) upwx
00135 real(8) upwy
00136 
00137 real(8) betaAC
00138 
00139 real(8) sourceT
00140 real(8) sourceU
00141 
00142 real(8) dx,dy
00143 
00144 real(8),pointer :: myLx,myLy
00145 
00146 
00147 real(8),pointer :: phyOrder
00148 
00149 
00150 real(8),pointer :: Re
00151 
00152 
00153 !< if x+eps*v is to be evaluated, then tempy is updated before calling discrete(1)
00154 !< such that tempy(setVec)=tempy(setVec)+lambda*deltax
00155 !< if discrete(0) is called, original y is evaluated
00156 do i=1,a
00157         discreteY=>tempy
00158 end do
00159 
00160 
00161 
00162 
00163 
00164 betaAC=10.d0
00165 
00166 !betaAC=1.d0/dTau
00167 
00168 discrete=0.d0
00169 
00170 
00171 x=>discreteY
00172 
00173 xPTO=>yPTO
00174 
00175 
00176 discrete=x(setVEc)
00177 
00178 
00179 
00180 ! In this problem, discretization for 1 and 2 is the same although they are located in different
00181 ! part of the domains of different physics
00182 ! they might be coded in a loop so double explicit definition for P1 and P2 is avoided
00183 ! this is not the case for BC 3 since it is defined explicitly
00184 ! however, if and interpolation interface is used, than the bourden is passed to the interface routine
00185 ! so that the discretizations looks the same. After then, P1 and P2 can be defined in a loop.
00186 
00187 
00188 
00189 
00190 
00191 
00192 
00193 
00194 !!!!!!!!!!!!!!
00195 !< assume that  physics-1 starts here
00196 !< also that the domain-1 of physics-1 starts here, as well
00197 !#############################################################################################
00198 !#############################################################################################
00199 !#############################################################################################
00200 !
00201 !############                                   ###
00202 !#############                                  ###
00203 !###            ###                                     ###
00204 !###            ###                                     ###
00205 !###            ###                                     ###
00206 !###            ###                                     ###
00207 !#############          ####            ###
00208 !############           ####            ###
00209 !###                                                    ###
00210 !###                                                    ###
00211 !###                                                    ###
00212 !###                                                    ###
00213 !###                                                    ###
00214 !###                                                    ###     
00215 !
00216 !#############################################################################################
00217 !#############################################################################################
00218 !#############################################################################################
00219 
00220 
00221 !modelPhysics=currentPhysics
00222 !modelDomain=currentDomain
00223 
00224 modelPhysics=1
00225 
00226 
00227 !< the idea is  mymodel%discrete%set(1)
00228 !<                              mymodel%discrete%set(2)
00229 !<                              mymodel%discrete%set(3)
00230 ! with set vec pointing to the set of the physics(..)%domain%(..)%set
00231 
00232 
00233 
00234 phyOrder=>myModel%physics(modelPhysics)%evalOrder
00235 
00236 RE=>myModel%physics(modelPhysics)%physicalParameters(1)
00237 
00238 
00239 mu=1.d0
00240 rho=1000.d0
00241 
00242 
00243 !Re=100.d0
00244 !Re=1000.d0
00245 !Re=2000.d0
00246 !Re=3000.d0
00247 !Re=3500.d0
00248 !Re=5000.d0
00249 !Re=10000.d0
00250 !stop
00251 
00252 
00253 
00254 !print*,'dx,dy=',dx,dy
00255 !stop
00256 
00257 
00258 dof=>myModel%physics(modelPhysics)%dof
00259 
00260 
00261 
00262 
00263 
00264 
00265 
00266 allocate(RHS(dof))
00267 RHS=0.d0
00268 
00269 !#####################################################################
00270 !user supplied variables
00271 !#####################################################################
00272 
00273 do myP=1,myModel%endIndex(modelPhysics)
00274 
00275         !#############################################################################################
00276         !
00277         !########                       ##
00278         !##             ##                      ##
00279         !##             ##                      ##
00280         !##             ##      ###             ##
00281         !##             ##                      ##
00282         !##             ##                      ##
00283         !########                       ##
00284         !
00285         !#############################################################################################
00286         
00287         modelDomain=1
00288 
00289         IM=>myModel%physics(modelPhysics)%domain(modelDomain)%dofMatrix
00290 
00291         myLx=>myModel%physics(modelPhysics)%domain(modelDomain)%Lx
00292         myLy=>myModel%physics(modelPhysics)%domain(modelDomain)%Ly
00293 
00294         !dx=myLx/(myModel%physics(modelPhysics)%domain(modelDomain)%maxI-myModel%physics(modelPhysics)%domain(modelDomain)%minI)
00295         !dy=myLy/(myModel%physics(modelPhysics)%domain(modelDomain)%maxJ-myModel%physics(modelPhysics)%domain(modelDomain)%minJ)
00296         dx=myLx/(myModel%physics(modelPhysics)%domain(modelDomain)%nx-1.d0)
00297         dy=myLy/(myModel%physics(modelPhysics)%domain(modelDomain)%ny-1.d0)
00298         ! normally, especially in finite volume coordinates should be read so that face lengths and center-to-center distances
00299         ! are calculated. for this uniform finite difference grid dx,dy is enough to proceed
00300 
00301         ! caution, how is dx dy in finite volume!!! was lx for the whole domain with ghostcells 
00302         ! or just the real interior domain w/o ghost cells!!! - if so number of ghostcells connection!!
00303 
00304         !print*,dx,dy
00305 
00306 
00307 
00308         do myD=1,myModel%physics(modelPhysics)%endIndex(modelDomain)
00309 
00310 
00311         !---------------------------------------------------------------------------------------------
00312         !interior
00313         !bc=1
00314         myBC=1
00315 
00316         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00317         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00318 
00319         do myJ=1,maxJ
00320                 
00321                 j=jVec(myJ)
00322                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00323                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00324 
00325                 RHS=0
00326                 do myI=1,maxI
00327 
00328                         i=iVec(myI)
00329                         
00330                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00331                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00332                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00333 
00334                         !unk1
00335                         t=1
00336                         do myT=1,endIndex(t)
00337 
00338                                 PC=x(IM(i,j,t,absolute))
00339                                 !PW=x(IM(i-1,j,t,absolute))
00340                                 !PE=x(IM(i+1,j,t,absolute))
00341                                 !PN=x(IM(i,j+1,t,absolute))
00342                                 !PS=x(IM(i,j-1,t,absolute))
00343 
00344 
00345                                 UC=x(IM(i,j,t+1,absolute))
00346                                 UW=x(IM(i-1,j,t+1,absolute))
00347                                 UE=x(IM(i+1,j,t+1,absolute))
00348                                 
00349                                 VC=x(IM(i,j,t+2,absolute))
00350                                 VS=x(IM(i,j-1,t+2,absolute))
00351                                 VN=x(IM(i,j+1,t+2,absolute))
00352 
00353                                 PCPTO=xPTO(IM(i,j,t,absolute))
00354 
00355 
00356                                 !---------
00357 
00358                                 dPdTau=(PC-PCPTO)*dx*dy/dTau
00359 
00360                                 !dPdX=(PE-PW)/(2.d0*dx)
00361                                 
00362                                 !dPdY=(PN-PS)/(2.d0*dy)
00363 
00364                                 rhoC=rho
00365                                 rhoW=rho
00366                                 rhoE=rho
00367                                 rhoN=rho
00368                                 rhoS=rho
00369                                 
00370 
00371                                 rhofw=(rhoc+rhow)/2.d0
00372                                 rhofe=(rhoc+rhoe)/2.d0
00373                                 rhofn=(rhoc+rhon)/2.d0
00374                                 rhofs=(rhoc+rhos)/2.d0
00375 
00376 
00377                                 Ufw=(Uc+Uw)/2.d0
00378                                 Ufe=(Uc+Ue)/2.d0
00379 
00380                                 Vfn=(Vc+Vn)/2.d0
00381                                 Vfs=(Vc+Vs)/2.d0
00382 
00383 
00384                                 !RHS(t)=dPdTau/betaAC+(Ufe*dy-Ufw*dy+Vfn*dx-Vfs*dx)
00385                                 RHS(t)=dPdTau/betaAC+(rhofe*Ufe*dy-rhofw*Ufw*dy+rhofn*Vfn*dx-rhofs*Vfs*dx)
00386 
00387                         end do
00388 
00389                         !unk2
00390                         t=2
00391                         do myT=1,endIndex(t)
00392 
00393                                 PW=x(IM(i-1,j,t-1,absolute))
00394                                 PE=x(IM(i+1,j,t-1,absolute))
00395 
00396                                 UC=x(IM(i,j,t,absolute))
00397                                 UW=x(IM(i-1,j,t,absolute))
00398                                 UE=x(IM(i+1,j,t,absolute))
00399                                 US=x(IM(i,j-1,t,absolute))
00400                                 UN=x(IM(i,j+1,t,absolute))
00401                                 UEE=x(IM(i+2,j,t,absolute))
00402                                 UWW=x(IM(i-2,j,t,absolute))
00403                                 UNN=x(IM(i,j+2,t,absolute))
00404                                 USS=x(IM(i,j-2,t,absolute))
00405                                 
00406                                 VC=x(IM(i,j,t+1,absolute))
00407                                 VN=x(IM(i,j+1,t+1,absolute))
00408                                 VS=x(IM(i,j-1,t+1,absolute))
00409 
00410                                 UCPTO=xPTO(IM(i,j,t,absolute))
00411 
00412                                 PCPTO=xPTO(IM(i,j,t-1,absolute))
00413 
00414 
00415                                 !---------
00416 
00417                                 dPdTau=(PC-PCPTO)*dx*dy/dTau
00418                                 
00419                                 !pseudo time
00420                                 dUdTau=(UC-UCPTO)*dx*dy/dTau
00421                                 
00422                                 !pressure gradient
00423                                 Px=(PE-PW)*dy/2.d0
00424 
00425                                 !diffusion
00426                                 dUdXw=(Uc-Uw)/dx
00427                                 dUdXe=(Ue-Uc)/dx
00428                                 dUdYs=(Uc-Us)/dy
00429                                 dUdYn=(Un-Uc)/dy
00430 
00431                                 U_fd_e=dy*dUdXe
00432                                 U_fd_w=dy*dUdXw
00433                                 U_fd_n=dx*dUdYn
00434                                 U_fd_s=dx*dUdYs
00435 
00436                                 U_fd=U_fd_e-U_fd_w+U_fd_n-U_fd_s
00437 
00438 
00439                                 !convection
00440 
00441                                 Ufw=(Uc+Uw)/2.d0
00442                                 Ufe=(Uc+Ue)/2.d0
00443                                 Ufn=(Uc+Un)/2.d0
00444                                 Ufs=(Uc+Us)/2.d0
00445 
00446                                 Vfn=(Vc+Vn)/2.d0
00447                                 Vfs=(Vc+Vs)/2.d0
00448 
00449 
00450 
00451                                 rhoC=rho
00452                                 rhoW=rho
00453                                 rhoE=rho
00454                                 rhoN=rho
00455                                 rhoS=rho
00456                                 
00457 
00458                                 rhofw=(rhoc+rhow)/2.d0
00459                                 rhofe=(rhoc+rhoe)/2.d0
00460                                 rhofn=(rhoc+rhon)/2.d0
00461                                 rhofs=(rhoc+rhos)/2.d0
00462 
00463 
00464 
00465 
00466 
00467                                 !-----------------------------------------------
00468                                 !###############################################
00469                                 !-----------------------------------------------
00470                                 !U_fc_e
00471 
00472                                 if (Ufe.gt.0) then 
00473 
00474                                         down=Ue
00475                                         upw1=Uc
00476                                         upw2=Uw
00477 
00478                                 else
00479 
00480                                         down=Uc
00481                                         upw1=Ue
00482                                         upw2=UEE
00483                                 
00484                                 end if
00485                                                                 
00486                                 value=(3.d0*down+6.d0*upw1-1.d0*upw2)/8.d0
00487 
00488                                 U_fc_e=value*Ufe
00489                                 
00490                                 U_fc_e=phyOrder*U_fc_e+(1.d0-phyOrder)*Ufe*Ufe
00491                                 !-----------------------------------------------
00492                                 !###############################################
00493                                 !-----------------------------------------------
00494                                 !U_fc_w
00495 
00496 
00497                                 if (Ufw.gt.0) then 
00498 
00499                                         down=Uc
00500                                         upw1=Uw
00501                                         upw2=UWW
00502 
00503                                 else
00504 
00505                                         down=Uw
00506                                         upw1=Uc
00507                                         upw2=Ue
00508 
00509                                 end if
00510 
00511                                 value=(3.d0*down+6.d0*upw1-1.d0*upw2)/8.d0
00512                                 
00513                                 U_fc_w=value*Ufw
00514 
00515                                 U_fc_w=phyOrder*U_fc_w+(1.d0-phyOrder)*Ufw*Ufw
00516                                 !-----------------------------------------------
00517                                 !###############################################
00518                                 !-----------------------------------------------
00519                                 !U_fc_n
00520                                 
00521 
00522                                 if (Vfn.gt.0) then 
00523 
00524                                         down=Un
00525                                         upw1=Uc
00526                                         upw2=Us
00527 
00528                                 else
00529 
00530                                         down=Uc
00531                                         upw1=Un
00532                                         upw2=UNN
00533 
00534                                 end if
00535                                 
00536 
00537                                 value=(3.d0*down+6.d0*upw1-1.d0*upw2)/8.d0
00538 
00539                                 U_fc_n=value*Vfn
00540 
00541                                 U_fc_n=phyOrder*U_fc_n+(1.d0-phyOrder)*Vfn*Ufn
00542                                 !-----------------------------------------------
00543                                 !###############################################
00544                                 !-----------------------------------------------
00545                                 !U_fc_s
00546 
00547 
00548                                 if (Vfs.gt.0) then 
00549 
00550                                         down=Uc
00551                                         upw1=Us
00552                                         upw2=USS
00553 
00554                                 else
00555 
00556                                         down=Us
00557                                         upw1=Uc
00558                                         upw2=Un
00559 
00560                                 end if
00561 
00562 
00563                                 value=(3.d0*down+6.d0*upw1-1.d0*upw2)/8.d0
00564 
00565                                 U_fc_s=value*Vfs
00566 
00567                                 U_fc_s=phyOrder*U_fc_s+(1.d0-phyOrder)*Vfs*Ufs
00568 
00569                                 !-----------------------------------------------
00570                                 !###############################################
00571                                 !-----------------------------------------------
00572 
00573                         
00574                                 !U_fc=U_fc_e*dy-U_fc_w*dy+U_fc_n*dx-U_fc_s*dx
00575                                 U_fc=rhofe*U_fc_e*dy-rhofw*U_fc_w*dy+rhofn*U_fc_n*dx-rhofs*U_fc_s*dx
00576 
00577 
00578                                 au=0.d0
00579                                 
00580                                 !betaAC=(max(dsqrt(10.d0),5.d0*(Uc*Uc+Vc*Vc)))**2.d0
00581 
00582                                 muC=mu
00583 
00584                                 !RHS(t)=(au*UC*dPdTau/betaAC+dUdTau+U_fc+Px)*Re-U_fd
00585                                 RHS(t)=(au*UC*dPdTau/betaAC+rhoC*dUdTau+U_fc+Px)-U_fd*muC
00586 
00587 
00588                         end do
00589 
00590                         !unk1
00591                         t=3
00592                         do myT=1,endIndex(t)
00593 
00594                                 PC=x(IM(i,j,t-2,absolute))
00595                                 PN=x(IM(i,j+1,t-2,absolute))
00596                                 PS=x(IM(i,j-1,t-2,absolute))
00597 
00598                                 UC=x(IM(i,j,t-1,absolute))
00599                                 UW=x(IM(i-1,j,t-1,absolute))
00600                                 UE=x(IM(i+1,j,t-1,absolute))
00601 
00602                                 VC=x(IM(i,j,t,absolute))
00603                                 VW=x(IM(i-1,j,t,absolute))
00604                                 VE=x(IM(i+1,j,t,absolute))
00605                                 VS=x(IM(i,j-1,t,absolute))
00606                                 VN=x(IM(i,j+1,t,absolute))
00607                                 VWW=x(IM(i-2,j,t,absolute))
00608                                 VEE=x(IM(i+2,j,t,absolute))
00609                                 VSS=x(IM(i,j-2,t,absolute))
00610                                 VNN=x(IM(i,j+2,t,absolute))
00611                                 
00612                                 VCPTO=xPTO(IM(i,j,t,absolute))
00613 
00614                                 PCPTO=xPTO(IM(i,j,t-2,absolute))
00615 
00616 
00617                                 !---------
00618 
00619                                 dPdTau=(PC-PCPTO)*dx*dy/dTau
00620                                 !pseudo time
00621                                 dVdTau=(VC-VCPTO)*dx*dy/dTau
00622 
00623                                 !pressure gradient
00624                                 Py=(PN-PS)*dx/2.d0
00625 
00626                                 !diffusion
00627                                 dVdXw=(Vc-Vw)/dx
00628                                 dVdXe=(Ve-Vc)/dx
00629                                 dVdYs=(Vc-Vs)/dy
00630                                 dVdYn=(Vn-Vc)/dy
00631 
00632                                 V_fd_e=dy*dVdXe
00633                                 V_fd_w=dy*dVdXw
00634                                 V_fd_n=dx*dVdYn
00635                                 V_fd_s=dx*dVdYs
00636 
00637                                 V_fd=V_fd_e-V_fd_w+V_fd_n-V_fd_s
00638 
00639                                 !convection
00640 
00641                                 Ufw=(Uc+Uw)/2.d0
00642                                 Ufe=(Uc+Ue)/2.d0
00643 
00644                                 Vfw=(Vc+Vw)/2.d0
00645                                 Vfe=(Vc+Ve)/2.d0
00646                                 Vfn=(Vc+Vn)/2.d0
00647                                 Vfs=(Vc+Vs)/2.d0
00648 
00649                                 rhoC=rho
00650                                 rhoW=rho
00651                                 rhoE=rho
00652                                 rhoN=rho
00653                                 rhoS=rho
00654                                 
00655 
00656                                 rhofw=(rhoc+rhow)/2.d0
00657                                 rhofe=(rhoc+rhoe)/2.d0
00658                                 rhofn=(rhoc+rhon)/2.d0
00659                                 rhofs=(rhoc+rhos)/2.d0
00660 
00661                                 !-----------------------------------------------
00662                                 !###############################################
00663                                 !-----------------------------------------------
00664                                 !V_fc_e
00665 
00666                                 if (Ufe.gt.0) then 
00667 
00668                                         down=Ve
00669                                         upw1=Vc
00670                                         upw2=Vw
00671 
00672 
00673                                 else
00674 
00675                                         down=Vc
00676                                         upw1=Ve
00677                                         upw2=VEE
00678 
00679                                 end if
00680 
00681 
00682                                 value=(3.d0*down+6.d0*upw1-1.d0*upw2)/8.d0
00683 
00684                                 V_fc_e=value*Ufe
00685 
00686                                 V_fc_e=phyOrder*V_fc_e+(1.d0-phyOrder)*Ufe*Vfe
00687 
00688 
00689                                 !-----------------------------------------------
00690                                 !###############################################
00691                                 !-----------------------------------------------
00692                                 !V_fc_w
00693 
00694                                 if (Ufw.gt.0) then 
00695 
00696                                         down=Vc
00697                                         upw1=Vw
00698                                         upw2=VWW
00699 
00700                                 else
00701 
00702                                         down=Vw
00703                                         upw1=Vc
00704                                         upw2=Ve
00705 
00706                                 end if
00707 
00708                                 value=(3.d0*down+6.d0*upw1-1.d0*upw2)/8.d0
00709 
00710                                 V_fc_w=value*Ufw
00711 
00712                                 V_fc_w=phyOrder*V_fc_w+(1.d0-phyOrder)*Ufw*Vfw
00713 
00714                                 !-----------------------------------------------
00715                                 !###############################################
00716                                 !-----------------------------------------------
00717                                 !V_fc_n
00718 
00719                                 if (Vfn.gt.0) then 
00720 
00721                                         down=Vn
00722                                         upw1=Vc
00723                                         upw2=Vs
00724 
00725                                 else
00726 
00727                                         down=Vc
00728                                         upw1=Vn
00729                                         upw2=VNN
00730 
00731                                 end if
00732 
00733 
00734                                 value=(3.d0*down+6.d0*upw1-1.d0*upw2)/8.d0
00735 
00736                                 V_fc_n=value*Vfn
00737 
00738                                 V_fc_n=phyOrder*V_fc_n+(1.d0-phyOrder)*Vfn*Vfn
00739 
00740                                 !-----------------------------------------------
00741                                 !###############################################
00742                                 !-----------------------------------------------
00743                                 !V_fc_s
00744 
00745                                 if (Vfs.gt.0) then 
00746 
00747                                         down=Vc
00748                                         upw1=Vs
00749                                         upw2=VSS
00750                                 else
00751 
00752                                         down=Vs
00753                                         upw1=Vc
00754                                         upw2=Vn
00755 
00756 
00757                                 end if
00758 
00759 
00760                                 value=(3.d0*down+6.d0*upw1-1.d0*upw2)/8.d0
00761 
00762                                 V_fc_s=value*Vfs
00763 
00764                                 V_fc_s=phyOrder*V_fc_s+(1.d0-phyOrder)*Vfs*Vfs
00765 
00766 
00767 
00768                                 !V_fc=V_fc_e*dy-V_fc_w*dy+V_fc_n*dx-V_fc_s*dx
00769                                 V_fc=rhofe*V_fc_e*dy-rhofw*V_fc_w*dy+rhofn*V_fc_n*dx-rhofs*V_fc_s*dx
00770 
00771 
00772                                 au=0.d0
00773 
00774                                 !betaAC=(max(dsqrt(10.d0),5.d0*(Uc*Uc+Vc*Vc)))**2.d0
00775 
00776                                 muC=mu
00777 
00778                                 !RHS(t)=(au*VC*dPdTau/betaAC+dVdTau+V_fc+Py)*Re-V_fd
00779                                 RHS(t)=(au*VC*dPdTau/betaAC+rhoC*dVdTau+V_fc+Py)-V_fd*muC
00780 
00781 
00782                         end do
00783 
00784 
00785 
00786                         do myT=1,maxT
00787 
00788                                 t=tVec(myT)
00789 
00790                                 !print*,IM(i,j,t,relative),rhs(t)
00791 
00792                                 discrete(IM(i,j,t,relative))=rhs(t)
00793 
00794                         end do
00795 
00796                 end do
00797 
00798 
00799         end do
00800         !---------------------------------------------------------------------------------------------
00801         !boundary
00802         !---------------------------------------------------------------------------------------------
00803         !bottom
00804         !bc=2
00805         myBC=2 
00806 
00807         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00808         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00809 
00810 
00811         do myJ=1,maxJ
00812                 
00813                 j=jVec(myJ)
00814                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00815                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00816 
00817                 RHS=0
00818                 do myI=1,maxI
00819 
00820                         i=iVec(myI)
00821                         
00822                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00823                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00824                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00825 
00826                         !unk1
00827                         t=1
00828                         do myT=1,endIndex(t)
00829 
00830                                 ! center values
00831                                 PC=x(IM(i,j,t,absolute))
00832                                 PN=x(IM(i,j+1,t,absolute))
00833 
00834                                 RHS(t)=(PN-PC)/dy
00835 
00836                         end do
00837 
00838                         !unk2
00839                         t=2
00840                         do myT=1,endIndex(t)
00841 
00842                                 ! center values
00843                                 UC=x(IM(i,j,t,absolute))
00844                                 UN=x(IM(i,j+1,t,absolute))
00845 
00846                                 !print*,TC,i,j
00847 
00848                                 RHS(t)=(UN+UC)/2.d0-0.d0
00849 
00850                         end do
00851 
00852                         !unk3
00853                         t=3
00854                         do myT=1,endIndex(t)
00855 
00856                                 ! center values
00857                                 VC=x(IM(i,j,t,absolute))
00858                                 VN=x(IM(i,j+1,t,absolute))
00859 
00860                                 !print*,TC,i,j
00861 
00862                                 RHS(t)=(VN+VC)/2.d0-0.d0
00863 
00864                         end do
00865 
00866 
00867                         do myT=1,maxT
00868                                 
00869                                 t=tVec(myT)
00870 
00871                                 !print*,IM(i,j,t,relative)
00872 
00873                                 discrete(IM(i,j,t,relative))=rhs(t)
00874 
00875 
00876                         end do
00877 
00878                 end do
00879 
00880 
00881         end do
00882         !---------------------------------------------------------------------------------------------
00883         !stop
00884         !---------------------------------------------------------------------------------------------
00885         !top
00886         !bc=3
00887         myBC=3 
00888 
00889         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00890         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00891 
00892 
00893         do myJ=1,maxJ
00894                 
00895                 j=jVec(myJ)
00896                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00897                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00898 
00899                 RHS=0
00900                 do myI=1,maxI
00901 
00902                         i=iVec(myI)
00903                         
00904                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00905                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00906                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00907 
00908                         !unk1
00909                         t=1
00910                         do myT=1,endIndex(t)
00911 
00912                                 ! center values
00913                                 PC=x(IM(i,j,t,absolute))
00914                                 PS=x(IM(i,j-1,t,absolute))
00915 
00916                                 RHS(t)=(PC-PS)/dy
00917 
00918                         end do
00919 
00920 
00921                         !unk2
00922                         t=2
00923                         do myT=1,endIndex(t)
00924 
00925                                 ! center values
00926                                 UC=x(IM(i,j,t,absolute))
00927                                 US=x(IM(i,j-1,t,absolute))
00928 
00929                                 !print*,TC,i,j
00930 
00931                                 RHS(t)=(UC+US)/2.d0-1.d0
00932 
00933                         end do
00934 
00935                         !unk3
00936                         t=3
00937                         do myT=1,endIndex(t)
00938 
00939                                 ! center values
00940                                 VC=x(IM(i,j,t,absolute))
00941                                 VS=x(IM(i,j-1,t,absolute))
00942 
00943                                 !print*,TC,i,j
00944 
00945                                 RHS(t)=(VC+VS)/2.d0-0.d0
00946 
00947                         end do
00948 
00949 
00950                         do myT=1,maxT
00951                                 
00952                                 t=tVec(myT)
00953 
00954                                 !print*,IM(i,j,t,relative)
00955 
00956                                 discrete(IM(i,j,t,relative))=rhs(t)
00957 
00958 
00959                         end do
00960 
00961                 end do
00962 
00963 
00964         end do
00965         !---------------------------------------------------------------------------------------------
00966         !stop
00967         !---------------------------------------------------------------------------------------------
00968         !left
00969         !bc=4
00970         myBC=4   
00971 
00972         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00973         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00974 
00975 
00976         do myJ=1,maxJ
00977                 
00978                 j=jVec(myJ)
00979                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00980                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00981 
00982                 RHS=0
00983                 do myI=1,maxI
00984 
00985                         i=iVec(myI)
00986                         
00987                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00988                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00989                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00990 
00991 
00992 
00993                         !unk1
00994                         t=1
00995                         do myT=1,endIndex(t)
00996 
00997                                 ! center values
00998                                 PC=x(IM(i,j,t,absolute))
00999                                 PE=x(IM(i+1,j,t,absolute))
01000                                 
01001                                 RHS(t)=(PE-PC)/dx
01002 
01003                         end do
01004 
01005                         !unk2
01006                         t=2
01007                         do myT=1,endIndex(t)
01008 
01009                                 ! center values
01010                                 UC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
01011                                 UE=x(IM(i+1,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
01012 
01013                                 RHS(t)=(UE+UC)/2.d0-0.d0
01014 
01015                         end do
01016 
01017                         !unk3
01018                         t=3
01019                         do myT=1,endIndex(t)
01020 
01021                                 ! center values
01022                                 VC=x(IM(i,j,t,absolute))
01023                                 VE=x(IM(i+1,j,t,absolute))
01024 
01025                                 RHS(t)=(VE+VC)/2.d0-0.d0
01026 
01027                         end do
01028 
01029 
01030                         do myT=1,maxT
01031                                 
01032                                 t=tVec(myT)
01033 
01034                                 !print*,IM(i,j,t,relative)
01035 
01036                                 discrete(IM(i,j,t,relative))=rhs(t)
01037 
01038                         end do
01039 
01040                 end do
01041 
01042 
01043 
01044         end do
01045         !---------------------------------------------------------------------------------------------
01046         !stop
01047         !---------------------------------------------------------------------------------------------
01048         !right
01049         !bc=5
01050         myBC=5 
01051 
01052         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01053         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01054 
01055 
01056         do myJ=1,maxJ
01057                 
01058                 j=jVec(myJ)
01059                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01060                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01061 
01062                 RHS=0
01063                 do myI=1,maxI
01064 
01065                         i=iVec(myI)
01066                         
01067                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01068                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01069                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01070 
01071                         !unk1
01072                         t=1
01073                         do myT=1,endIndex(t)
01074 
01075                                 ! center values
01076                                 PC=x(IM(i,j,t,absolute))
01077                                 PW=x(IM(i-1,j,t,absolute))
01078 
01079                                 !print*,TC,i,j
01080 
01081                                 RHS(t)=(PC-PW)/dx
01082 
01083                         end do
01084 
01085                         !unk2
01086                         t=2
01087                         do myT=1,endIndex(t)
01088 
01089                                 ! center values
01090                                 UC=x(IM(i,j,t,absolute))
01091                                 UW=x(IM(i-1,j,t,absolute))
01092 
01093                                 !print*,TC,i,j
01094 
01095                                 RHS(t)=(UC+UW)/2.d0-0.d0
01096 
01097                         end do
01098 
01099                         !unk3
01100                         t=3
01101                         do myT=1,endIndex(t)
01102 
01103                                 ! center values
01104                                 VC=x(IM(i,j,t,absolute))
01105                                 VW=x(IM(i-1,j,t,absolute))
01106 
01107                                 !print*,TC,i,j
01108 
01109                                 RHS(t)=(VC+VW)/2.d0-0.d0
01110 
01111                         end do
01112 
01113 
01114                         do myT=1,maxT
01115                                 
01116                                 t=tVec(myT)
01117 
01118                                 !print*,IM(i,j,t,relative)
01119 
01120                                 discrete(IM(i,j,t,relative))=rhs(t)
01121 
01122 
01123                         end do
01124 
01125                 end do
01126 
01127 
01128         end do
01129         !---------------------------------------------------------------------------------------------
01130         !stop
01131         !---------------------------------------------------------------------------------------------
01132         !bottom - ghost cell 2
01133         !bc=6
01134         myBC=6 
01135 
01136         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01137         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01138 
01139 
01140         do myJ=1,maxJ
01141                 
01142                 j=jVec(myJ)
01143                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01144                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01145 
01146                 RHS=0
01147                 do myI=1,maxI
01148 
01149                         i=iVec(myI)
01150                         
01151                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01152                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01153                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01154 
01155                         !unk1
01156                         t=1
01157                         do myT=1,endIndex(t)
01158 
01159                                 ! center values
01160                                 PC=x(IM(i,j,t,absolute))
01161                                 PN=x(IM(i,j+1,t,absolute))
01162 
01163                                 RHS(t)=(PN-PC)/dy
01164 
01165                         end do
01166 
01167                         !unk2
01168                         t=2
01169                         do myT=1,endIndex(t)
01170 
01171                                 ! center values
01172                                 UC=x(IM(i,j,t,absolute))
01173                                 UN=x(IM(i,j+1,t,absolute))
01174                                 UNNN=x(IM(i,j+3,t,absolute))
01175 
01176                                 RHS(t)=phyOrder*((UNNN+UC)/2.d0-0.d0) + (1.d0-phyOrder)*(UC-3.d0*UN+2.d0*0.d0)
01177 
01178                         end do
01179 
01180                         !unk3
01181                         t=3
01182                         do myT=1,endIndex(t)
01183 
01184                                 ! center values
01185                                 VC=x(IM(i,j,t,absolute))
01186                                 VN=x(IM(i,j+1,t,absolute))
01187                                 VNNN=x(IM(i,j+3,t,absolute))
01188 
01189                                 RHS(t)=phyOrder*((VNNN+VC)/2.d0-0.d0) + (1.d0-phyOrder)*(VC-3.d0*VN+2.d0*0.d0)
01190 
01191                         end do
01192 
01193 
01194                         do myT=1,maxT
01195                                 
01196                                 t=tVec(myT)
01197 
01198                                 !print*,IM(i,j,t,relative)
01199 
01200                                 discrete(IM(i,j,t,relative))=rhs(t)
01201 
01202 
01203                         end do
01204 
01205                 end do
01206 
01207 
01208         end do
01209         !---------------------------------------------------------------------------------------------
01210         !stop
01211         !---------------------------------------------------------------------------------------------
01212         !top - ghost cell 2
01213         !bc=7
01214         myBC=7 
01215 
01216         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01217         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01218 
01219 
01220         do myJ=1,maxJ
01221                 
01222                 j=jVec(myJ)
01223                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01224                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01225 
01226                 RHS=0
01227                 do myI=1,maxI
01228 
01229                         i=iVec(myI)
01230                         
01231                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01232                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01233                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01234 
01235                         !unk1
01236                         t=1
01237                         do myT=1,endIndex(t)
01238 
01239                                 ! center values
01240                                 PC=x(IM(i,j,t,absolute))
01241                                 PS=x(IM(i,j-1,t,absolute))
01242 
01243                                 RHS(t)=(PC-PS)/dy
01244 
01245                         end do
01246 
01247 
01248                         !unk2
01249                         t=2
01250                         do myT=1,endIndex(t)
01251 
01252                                 ! center values
01253                                 UC=x(IM(i,j,t,absolute))
01254                                 US=x(IM(i,j-1,t,absolute))
01255                                 USSS=x(IM(i,j-3,t,absolute))
01256 
01257                                 RHS(t)=phyOrder*((USSS+UC)/2.d0-1.d0) + (1.d0-phyOrder)*(UC-3.d0*US+2.d0*1.d0)
01258 
01259                         end do
01260 
01261                         !unk3
01262                         t=3
01263                         do myT=1,endIndex(t)
01264 
01265                                 ! center values
01266                                 VC=x(IM(i,j,t,absolute))
01267                                 VS=x(IM(i,j-1,t,absolute))
01268                                 VSSS=x(IM(i,j-3,t,absolute))
01269 
01270                                 !print*,TC,i,j
01271 
01272                                 RHS(t)=phyOrder*((VSSS+VC)/2.d0-0.d0) + (1.d0-phyOrder)*(VC-3.d0*VS+2.d0*0.d0)
01273 
01274                         end do
01275 
01276 
01277                         do myT=1,maxT
01278                                 
01279                                 t=tVec(myT)
01280 
01281                                 !print*,IM(i,j,t,relative)
01282 
01283                                 discrete(IM(i,j,t,relative))=rhs(t)
01284 
01285 
01286                         end do
01287 
01288                 end do
01289 
01290 
01291         end do
01292         !---------------------------------------------------------------------------------------------
01293         !stop
01294         !---------------------------------------------------------------------------------------------
01295         !left - ghost cell 2
01296         !bc=8
01297         myBC=8   
01298 
01299         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01300         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01301 
01302 
01303         do myJ=1,maxJ
01304                 
01305                 j=jVec(myJ)
01306                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01307                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01308 
01309                 RHS=0
01310                 do myI=1,maxI
01311 
01312                         i=iVec(myI)
01313                         
01314                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01315                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01316                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01317 
01318 
01319 
01320                         !unk1
01321                         t=1
01322                         do myT=1,endIndex(t)
01323 
01324                                 ! center values
01325                                 PC=x(IM(i,j,t,absolute))
01326                                 PE=x(IM(i+1,j,t,absolute))
01327                                 
01328                                 RHS(t)=(PE-PC)/dx
01329 
01330                         end do
01331 
01332                         !unk2
01333                         t=2
01334                         do myT=1,endIndex(t)
01335 
01336                                 ! center values
01337                                 UC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
01338                                 UE=x(IM(i+1,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
01339                                 UEEE=x(IM(i+3,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
01340 
01341                                 RHS(t)=phyOrder*((UEEE+UC)/2.d0-0.d0) + (1.d0-phyOrder)*(UC-3.d0*UE+2.d0*0.d0)
01342 
01343                         end do
01344 
01345                         !unk3
01346                         t=3
01347                         do myT=1,endIndex(t)
01348 
01349                                 ! center values
01350                                 VC=x(IM(i,j,t,absolute))
01351                                 VE=x(IM(i+1,j,t,absolute))
01352                                 VEEE=x(IM(i+3,j,t,absolute))
01353 
01354                                 RHS(t)=phyOrder*((VEEE+VC)/2.d0-0.d0) + (1.d0-phyOrder)*(VC-3.d0*VE+2.d0*0.d0)
01355 
01356                         end do
01357 
01358 
01359                         do myT=1,maxT
01360                                 
01361                                 t=tVec(myT)
01362 
01363                                 !print*,IM(i,j,t,relative)
01364 
01365                                 discrete(IM(i,j,t,relative))=rhs(t)
01366 
01367                         end do
01368 
01369                 end do
01370 
01371 
01372 
01373         end do
01374         !---------------------------------------------------------------------------------------------
01375         !stop
01376         !---------------------------------------------------------------------------------------------
01377         !right - ghost cell 2
01378         !bc=9
01379         myBC=9 
01380 
01381         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
01382         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
01383 
01384 
01385         do myJ=1,maxJ
01386                 
01387                 j=jVec(myJ)
01388                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
01389                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
01390 
01391                 RHS=0
01392                 do myI=1,maxI
01393 
01394                         i=iVec(myI)
01395                         
01396                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
01397                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
01398                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
01399 
01400                         !unk1
01401                         t=1
01402                         do myT=1,endIndex(t)
01403 
01404                                 ! center values
01405                                 PC=x(IM(i,j,t,absolute))
01406                                 PW=x(IM(i-1,j,t,absolute))
01407 
01408                                 !print*,TC,i,j
01409 
01410                                 RHS(t)=(PC-PW)/dx
01411 
01412                         end do
01413 
01414                         !unk2
01415                         t=2
01416                         do myT=1,endIndex(t)
01417 
01418                                 ! center values
01419                                 UC=x(IM(i,j,t,absolute))
01420                                 UW=x(IM(i-1,j,t,absolute))
01421                                 UWWW=x(IM(i-3,j,t,absolute))
01422 
01423                                 !print*,TC,i,j
01424 
01425                                 RHS(t)=phyOrder*((UWWW+UC)/2.d0-0.d0) + (1.d0-phyOrder)*(UC-3.d0*UW+2.d0*0.d0)
01426 
01427                         end do
01428 
01429                         !unk3
01430                         t=3
01431                         do myT=1,endIndex(t)
01432 
01433                                 ! center values
01434                                 VC=x(IM(i,j,t,absolute))
01435                                 VW=x(IM(i-1,j,t,absolute))
01436                                 VWWW=x(IM(i-3,j,t,absolute))
01437 
01438                                 !print*,TC,i,j
01439 
01440                                 RHS(t)=phyOrder*((VWWW+VC)/2.d0-0.d0) + (1.d0-phyOrder)*(VC-3.d0*VW+2.d0*0.d0)
01441 
01442                         end do
01443 
01444 
01445                         do myT=1,maxT
01446                                 
01447                                 t=tVec(myT)
01448 
01449                                 !print*,IM(i,j,t,relative)
01450 
01451                                 discrete(IM(i,j,t,relative))=rhs(t)
01452 
01453 
01454                         end do
01455 
01456                 end do
01457 
01458 
01459         end do
01460         !---------------------------------------------------------------------------------------------
01461 
01462         end do ! domain 1 loop
01463 
01464 
01465 
01466 
01467 end do ! physics loop
01468 
01469 !#####################################################################
01470 
01471 
01472 deallocate(RHS)
01473 
01474 
01475 
01476 
01477 
01478 
01479 !===================================================================
01480 discreteY=>y
01481 
01482 end function discrete
01483 
01484 
01485 
01486 
01487 
01488 
01489 
01490 
01491 
01492 
01493 
01494 
01495 end module model
D:/demona/2.0/manual/doc/2-31set/problems/ldc_fv/1/modelMLN1Dim.f90
