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) VNN,VSS
00061 real(8) VNNN,VSSS
00062 
00063 
00064 real(8) PCPTO
00065 
00066 real(8) UCPTO,UCRTO
00067 
00068 real(8) VCPTO,VCRTO
00069 
00070 
00071 real(8) dPdX,dPdY
00072 real(8) dUdX,dUdY
00073 real(8) dVdX,dVdY
00074 real(8) d2UdX2,d2UdY2
00075 real(8) d2VdX2,d2VdY2
00076 
00077 real(8) dPdTau
00078 
00079 real(8) dUdTau,dUdT
00080 
00081 real(8) dVdTau,dVdT
00082 
00083 
00084 
00085 real(8) UdUdX,VdUdY
00086 real(8) UdVdX,VdVdY
00087 
00088 
00089 !real(8) dtau
00090 
00091 real(8) upu,upv
00092 
00093 real(8) myU,myV
00094 
00095 real(8) upwx
00096 real(8) upwy
00097 
00098 real(8) betaAC
00099 
00100 real(8) sourceT
00101 real(8) sourceU
00102 
00103 real(8) dx,dy
00104 
00105 real(8),pointer :: myLx,myLy
00106 
00107 
00108 real(8),pointer :: phyOrder
00109 
00110 
00111 real(8),pointer :: Re
00112 
00113 
00114 !< if x+eps*v is to be evaluated, then tempy is updated before calling discrete(1)
00115 !< such that tempy(setVec)=tempy(setVec)+lambda*deltax
00116 !< if discrete(0) is called, original y is evaluated
00117 do i=1,a
00118         discreteY=>tempy
00119 end do
00120 
00121 
00122 
00123 
00124 
00125 betaAC=10.d0
00126 
00127 
00128 
00129 discrete=0.d0
00130 
00131 
00132 x=>discreteY
00133 
00134 xPTO=>yPTO
00135 
00136 
00137 discrete=x(setVEc)
00138 
00139 
00140 
00141 ! In this problem, discretization for 1 and 2 is the same although they are located in different
00142 ! part of the domains of different physics
00143 ! they might be coded in a loop so double explicit definition for P1 and P2 is avoided
00144 ! this is not the case for BC 3 since it is defined explicitly
00145 ! however, if and interpolation interface is used, than the bourden is passed to the interface routine
00146 ! so that the discretizations looks the same. After then, P1 and P2 can be defined in a loop.
00147 
00148 
00149 
00150 
00151 
00152 
00153 
00154 
00155 !!!!!!!!!!!!!!
00156 !< assume that  physics-1 starts here
00157 !< also that the domain-1 of physics-1 starts here, as well
00158 !#############################################################################################
00159 !#############################################################################################
00160 !#############################################################################################
00161 !
00162 !############                                   ###
00163 !#############                                  ###
00164 !###            ###                                     ###
00165 !###            ###                                     ###
00166 !###            ###                                     ###
00167 !###            ###                                     ###
00168 !#############          ####            ###
00169 !############           ####            ###
00170 !###                                                    ###
00171 !###                                                    ###
00172 !###                                                    ###
00173 !###                                                    ###
00174 !###                                                    ###
00175 !###                                                    ###     
00176 !
00177 !#############################################################################################
00178 !#############################################################################################
00179 !#############################################################################################
00180 
00181 
00182 !modelPhysics=currentPhysics
00183 !modelDomain=currentDomain
00184 
00185 modelPhysics=1
00186 
00187 
00188 !< the idea is  mymodel%discrete%set(1)
00189 !<                              mymodel%discrete%set(2)
00190 !<                              mymodel%discrete%set(3)
00191 ! with set vec pointing to the set of the physics(..)%domain%(..)%set
00192 
00193 
00194 
00195 phyOrder=>myModel%physics(modelPhysics)%evalOrder
00196 
00197 RE=>myModel%physics(modelPhysics)%physicalParameters(1)
00198 
00199 
00200 
00201 !Re=100.d0
00202 !Re=1000.d0
00203 !Re=2000.d0
00204 !Re=3000.d0
00205 !Re=3500.d0
00206 !Re=5000.d0
00207 !Re=10000.d0
00208 !stop
00209 
00210 
00211 
00212 !print*,'dx,dy=',dx,dy
00213 !stop
00214 
00215 
00216 dof=>myModel%physics(modelPhysics)%dof
00217 
00218 
00219 
00220 
00221 allocate(RHS(dof))
00222 RHS=0.d0
00223 
00224 !#####################################################################
00225 !user supplied variables
00226 !#####################################################################
00227 
00228 do myP=1,myModel%endIndex(modelPhysics)
00229 
00230         !#############################################################################################
00231         !
00232         !########                       ##
00233         !##             ##                      ##
00234         !##             ##                      ##
00235         !##             ##      ###             ##
00236         !##             ##                      ##
00237         !##             ##                      ##
00238         !########                       ##
00239         !
00240         !#############################################################################################
00241         
00242         modelDomain=1
00243 
00244         IM=>myModel%physics(modelPhysics)%domain(modelDomain)%dofMatrix
00245 
00246         myLx=>myModel%physics(modelPhysics)%domain(modelDomain)%Lx
00247         myLy=>myModel%physics(modelPhysics)%domain(modelDomain)%Ly
00248 
00249         dx=myLx/(myModel%physics(modelPhysics)%domain(modelDomain)%maxI-myModel%physics(modelPhysics)%domain(modelDomain)%minI)
00250         dy=myLy/(myModel%physics(modelPhysics)%domain(modelDomain)%maxJ-myModel%physics(modelPhysics)%domain(modelDomain)%minJ)
00251         ! normally, especially in finite volume coordinates should be read so that face lengths and center-to-center distances
00252         ! are calculated. for this uniform finite difference grid dx,dy is enough to proceed
00253 
00254         ! caution, how is dx dy in finite volume!!! was lx for the whole domain with ghostcells 
00255         ! or just the real interior domain w/o ghost cells!!! - if so number of ghostcells connection!!
00256 
00257         !print*,dx,dy
00258 
00259 
00260 
00261         do myD=1,myModel%physics(modelPhysics)%endIndex(modelDomain)
00262 
00263 
00264         !---------------------------------------------------------------------------------------------
00265         !interior
00266         !bc=1
00267         myBC=1
00268 
00269         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00270         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00271 
00272         do myJ=1,maxJ
00273                 
00274                 j=jVec(myJ)
00275                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00276                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00277 
00278                 RHS=0
00279                 do myI=1,maxI
00280 
00281                         i=iVec(myI)
00282                         
00283                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00284                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00285                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00286 
00287                         !unk1
00288                         t=1
00289                         do myT=1,endIndex(t)
00290 
00291                                 PC=x(IM(i,j,t,absolute))
00292                                 PW=x(IM(i-1,j,t,absolute))
00293                                 PE=x(IM(i+1,j,t,absolute))
00294                                 PN=x(IM(i,j+1,t,absolute))
00295                                 PS=x(IM(i,j-1,t,absolute))
00296 
00297 
00298                                 UC=x(IM(i,j,t+1,absolute))
00299                                 UW=x(IM(i-1,j,t+1,absolute))
00300                                 UE=x(IM(i+1,j,t+1,absolute))
00301                                 
00302                                 VC=x(IM(i,j,t+2,absolute))
00303                                 VS=x(IM(i,j-1,t+2,absolute))
00304                                 VN=x(IM(i,j+1,t+2,absolute))
00305 
00306                                 PCPTO=xPTO(IM(i,j,t,absolute))
00307 
00308                                 dPdTau=(PC-PCPTO)/(dTau)
00309 
00310                                 dPdX=(PE-PC)/(1.d0*dx)
00311                                 
00312                                 dPdY=(PN-PC)/(1.d0*dy)
00313 
00314                                 dUdX=(UE-UW)/(2.d0*dx)
00315 
00316                                 dVdY=(VN-VS)/(2.d0*dy)
00317 
00318 
00319                                 RHS(t)=dPdTau+0.d0*(UC*dPdX+VC*dPdY)+betaAC*(dUdX+dVdY)
00320 
00321                                 !print*,dPdTau +betaAC*(dUdX+dVdY)
00322                                 !stop
00323 
00324                                 !RHS(t)=Pc-i+j
00325 
00326                                 !RHS(t)=phyOrder*rhs(t)+(1.d0-phyOrder)*Pc/dTau  ! so we will have a diagonal preconditioner
00327 
00328 
00329                         end do
00330 
00331                         !unk2
00332                         t=2
00333                         do myT=1,endIndex(t)
00334 
00335 
00336                                 PC=x(IM(i,j,t-1,absolute))
00337                                 PW=x(IM(i-1,j,t-1,absolute))
00338                                 PE=x(IM(i+1,j,t-1,absolute))
00339 
00340                                 UC=x(IM(i,j,t,absolute))
00341                                 UW=x(IM(i-1,j,t,absolute))
00342                                 UE=x(IM(i+1,j,t,absolute))
00343                                 US=x(IM(i,j-1,t,absolute))
00344                                 UN=x(IM(i,j+1,t,absolute))
00345                                 
00346                                 VC=x(IM(i,j,t+1,absolute))
00347 
00348                                 UCPTO=xPTO(IM(i,j,t,absolute))
00349 
00350                                 dUdTau=(UC-UCPTO)/dTau
00351 
00352 
00353                                 dPdX=(PE-PW)/(2.d0*dx)
00354 
00355                                 dUdX=(UE-UW)/(2.d0*dx)
00356                                 dUdY=(UN-US)/(2.d0*dy)
00357 
00358 
00359                                 upwx=signum(Uc)
00360                                 upwy=signum(Vc)
00361 
00362 
00363                                 UdUdX=Uc*(1.d0-upwx)*(Ue-Uc)/(2.d0*dx)+Uc*(1.d0+upwx)*(Uc-Uw)/(2.d0*dx)
00364                                 VdUdY=Vc*(1.d0-upwy)*(Un-Uc)/(2.d0*dy)+Vc*(1.d0+upwy)*(Uc-Us)/(2.d0*dy)
00365 
00366                                 UdUdX=Uc*dUdX
00367                                 VdUdY=Vc*dUdY
00368 
00369                                 d2UdX2=(UW+UE-2.d0*UC)/(dx**2.d0)
00370                                 d2UdY2=(UN+US-2.d0*UC)/(dy**2.d0)
00371 
00372 
00373 
00374                                 RHS(t)=Re*(0.d0*uc/dt+dUdTau + UdUdX + VdUdY + dPdX) - d2UdX2 - d2UdY2
00375                                 !RHS(t)=- d2UdX2 - d2UdY2
00376 
00377                                 !RHS(t)=phyOrder*rhs(t)+(1.d0-phyOrder)*(2.d0*UC/(dx**2.d0)+2.d0*UC/(dy**2.d0))
00378 
00379 
00380                         end do
00381 
00382                         !unk1
00383                         t=3
00384                         do myT=1,endIndex(t)
00385 
00386                                 PC=x(IM(i,j,t-2,absolute))
00387                                 PN=x(IM(i,j+1,t-2,absolute))
00388                                 PS=x(IM(i,j-1,t-2,absolute))
00389 
00390                                 UC=x(IM(i,j,t-1,absolute))
00391 
00392                                 VC=x(IM(i,j,t,absolute))
00393                                 VW=x(IM(i-1,j,t,absolute))
00394                                 VE=x(IM(i+1,j,t,absolute))
00395                                 VS=x(IM(i,j-1,t,absolute))
00396                                 VN=x(IM(i,j+1,t,absolute))
00397                                 
00398                                 VCPTO=xPTO(IM(i,j,t,absolute))
00399 
00400                                 dVdTau=(VC-VCPTO)/dTau
00401 
00402 
00403                                 dPdY=(PN-PS)/(2.d0*dy)
00404 
00405                                 dVdX=(VE-VW)/(2.d0*dx)
00406                                 
00407                                 dVdY=(VN-VS)/(2.d0*dy)
00408 
00409 
00410 
00411                                 upwx=signum(Uc)
00412                                 upwy=signum(Vc)
00413 
00414                                 
00415                                 UdVdX=Uc*(1.d0-upwx)*(Ve-Vc)/(2.d0*dx)+Uc*(1.d0+upwx)*(Vc-Vw)/(2.d0*dx)
00416 
00417                                 VdVdY=Vc*(1.d0-upwy)*(Vn-Vc)/(2.d0*dy)+Vc*(1.d0+upwy)*(Vc-Vs)/(2.d0*dy)
00418 
00419                                 
00420                                 UdVdX=Uc*dVdX                   
00421                                 VdVdY=Vc*dVdY
00422 
00423 
00424 
00425                                 d2VdX2=(VW+VE-2.d0*VC)/(dx**2.d0)
00426 
00427                                 d2VdY2=(VN+VS-2.d0*VC)/(dy**2.d0)
00428 
00429 
00430                         
00431 
00432                                 RHS(t)=Re*(0.d0*vc/dt+dVdTau + UdVdX + VdVdY + dPdY) - d2VdX2 - d2VdY2
00433 
00434                                 !RHS(t)=Re*(dVdTau + UdVdX + VdVdY + dPdY) - d2VdX2 - d2VdY2+1.d0
00435 
00436 
00437                                 !RHS(t)=Re*(UC*dVORdX + VC*dVORdY)-(d2VORdX2+d2VORdY2)
00438                                 !RHS(t)=RE*(UdVORdX+VdVORdY)-(d2VORdX2+d2VORdY2)
00439 
00440                                 !RHS(t)=phyOrder*rhs(t)+(1.d0-phyOrder)*(2.d0*VC/(dx**2.d0)+2.d0*VC/(dy**2.d0))
00441 
00442 
00443                         end do
00444 
00445 
00446 
00447                         do myT=1,maxT
00448 
00449                                 t=tVec(myT)
00450 
00451                                 !print*,IM(i,j,t,relative),rhs(t)
00452 
00453                                 discrete(IM(i,j,t,relative))=rhs(t)
00454 
00455                         end do
00456 
00457                 end do
00458 
00459 
00460         end do
00461         !---------------------------------------------------------------------------------------------
00462         !boundary
00463         !---------------------------------------------------------------------------------------------
00464         !left
00465         !bc=2
00466         myBC=2   
00467 
00468         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00469         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00470 
00471 
00472         do myJ=1,maxJ
00473                 
00474                 j=jVec(myJ)
00475                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00476                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00477 
00478                 RHS=0
00479                 do myI=1,maxI
00480 
00481                         i=iVec(myI)
00482                         
00483                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00484                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00485                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00486 
00487 
00488 
00489                         !unk1
00490                         t=1
00491                         do myT=1,endIndex(t)
00492 
00493                                 ! center values
00494                                 PC=x(IM(i,j,t,absolute))
00495                                 PE=x(IM(i+1,j,t,absolute))
00496                                 PEE=x(IM(i+2,j,t,absolute))
00497 
00498                                 ! center values
00499                                 UC=x(IM(i,j,t+1,absolute))
00500                                 UE=x(IM(i+1,j,t+1,absolute))
00501                                 UEE=x(IM(i+2,j,t+1,absolute))
00502                                 UEEE=x(IM(i+3,j,t+1,absolute))
00503                                 
00504                                 !print*,TC,i,j
00505 
00506                                 !RHS(t)=VORC-(VE-VC)/dx
00507                                 !RHS(t)=PC-phyOrder*(-3.d0*VC+4.d0*VE-x(IM(i+2,j,t-1,absolute)))/(2.d0*dx)-(1.d0-phyOrder)*(VE-VC)/dx
00508                                 !RHS(t)=VORC-(-3.d0*VC+4.d0*VE-x(IM(i+2,j,t-1,absolute)))/(2.d0*dx)
00509 
00510 
00511 
00512                                 !RHS(t)=phyOrder*(Re*(-3.d0*PC+4.d0*PE-1.d0*PEE)/(2.d0*dx)-(2.d0*UC-5.d0*UE+4.d0*UEE-1.d0*UEEE)/dx**2.d0)+(1.d0-phyOrder)*(PE-PC)/dx
00513                                 !RHS(t)=phyOrder*(Re*(-3.d0*PC+4.d0*PE-1.d0*PEE)/(2.d0*dx)-(2.d0*UC-5.d0*UE+4.d0*UEE-1.d0*UEEE)/dx**2.d0)+(1.d0-phyOrder)*(Re*(-3.d0*PC+4.d0*PE)/(2.d0*dx)-(2.d0*UC-5.d0*UE)/dx**2.d0)
00514                                 RHS(t)=phyOrder*(Re*(-3.d0*PC+4.d0*PE-1.d0*PEE)/(2.d0*dx)-(2.d0*UC-5.d0*UE+4.d0*UEE-1.d0*UEEE)/dx**2.d0)+(1.d0-phyOrder)*(Re*(-3.d0*PC)/(2.d0*dx))
00515                                 RHS(t)=(Re*(-3.d0*PC+4.d0*PE-1.d0*PEE)/(2.d0*dx)-(2.d0*UC-5.d0*UE+4.d0*UEE-1.d0*UEEE)/dx**2.d0)
00516 
00517 
00518                         end do
00519 
00520                         !unk2
00521                         t=2
00522                         do myT=1,endIndex(t)
00523 
00524                                 ! center values
00525                                 UC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00526 
00527                                 !print*,TC,i,j
00528 
00529                                 RHS(t)=UC-0.d0
00530 
00531                         end do
00532 
00533                         !unk3
00534                         t=3
00535                         do myT=1,endIndex(t)
00536 
00537                                 ! center values
00538                                 VC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00539 
00540                                 !print*,TC,i,j
00541 
00542                                 RHS(t)=VC-0.d0
00543 
00544                         end do
00545 
00546 
00547                         do myT=1,maxT
00548                                 
00549                                 t=tVec(myT)
00550 
00551                                 !print*,IM(i,j,t,relative)
00552 
00553                                 discrete(IM(i,j,t,relative))=rhs(t)
00554 
00555                         end do
00556 
00557                 end do
00558 
00559 
00560 
00561         end do
00562         !---------------------------------------------------------------------------------------------
00563         !stop
00564         !---------------------------------------------------------------------------------------------
00565         !right
00566         !bc=3
00567         myBC=3 
00568 
00569         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00570         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00571 
00572 
00573         do myJ=1,maxJ
00574                 
00575                 j=jVec(myJ)
00576                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00577                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00578 
00579                 RHS=0
00580                 do myI=1,maxI
00581 
00582                         i=iVec(myI)
00583                         
00584                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00585                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00586                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00587 
00588                         !unk1
00589                         t=1
00590                         do myT=1,endIndex(t)
00591 
00592                                 ! center values
00593                                 PC=x(IM(i,j,t,absolute))
00594                                 PW=x(IM(i-1,j,t,absolute))
00595                                 PWW=x(IM(i-2,j,t,absolute))
00596                                 
00597                                 UC=x(IM(i,j,t+1,absolute))
00598                                 UW=x(IM(i-1,j,t+1,absolute))
00599                                 UWW=x(IM(i-2,j,t+1,absolute))
00600                                 UWWW=x(IM(i-3,j,t+1,absolute))
00601 
00602                                 !print*,TC,i,j
00603 
00604                                 !RHS(t)=VORC-(VC-VW)/dx
00605 
00606                                 !RHS(t)=PC-phyOrder*(3.d0*VC-4.d0*VW+x(IM(i-2,j,t-1,absolute)))/(2.d0*dx)-(1.d0-phyOrder)*(VC-VW)/dx
00607                                 !RHS(t)=VORC-(3.d0*VC-4.d0*VW+x(IM(i-2,j,t-1,absolute)))/(2.d0*dx)
00608 
00609                                 !RHS(t)=phyOrder*(Re*(3.d0*PC-4.d0*PW+1.d0*PWW)/(2.d0*dx)-(2.d0*UC-5.d0*UW+4.d0*UWW-1.d0*UWWW)/dx**2.d0) + (1.d0-phyOrder)*(PC-PW)/dx ! Re*del[p]/del[x]=del^2u/delx^2
00610                                 RHS(t)=phyOrder*(Re*(3.d0*PC-4.d0*PW+1.d0*PWW)/(2.d0*dx)-(2.d0*UC-5.d0*UW+4.d0*UWW-1.d0*UWWW)/dx**2.d0) + (1.d0-phyOrder)*(Re*(3.d0*PC-4.d0*PW)/(2.d0*dx)-(2.d0*UC-5.d0*UW)/dx**2.d0)
00611                                 RHS(t)=phyOrder*(Re*(3.d0*PC-4.d0*PW+1.d0*PWW)/(2.d0*dx)-(2.d0*UC-5.d0*UW+4.d0*UWW-1.d0*UWWW)/dx**2.d0) + (1.d0-phyOrder)*(Re*(3.d0*PC)/(2.d0*dx))
00612                                 RHS(t)=(Re*(3.d0*PC-4.d0*PW+1.d0*PWW)/(2.d0*dx)-(2.d0*UC-5.d0*UW+4.d0*UWW-1.d0*UWWW)/dx**2.d0)
00613 
00614                                 !rhs(t)=pc-0.d0
00615 
00616                                 if (test.eq.1) then
00617                                 
00618                                 !print*,PC,dx,phyOrder,(1.d0-phyOrder)*(3.d0*PC)/(2.d0*dx)
00619                                 
00620                                 !stop
00621 
00622                                 end if
00623 
00624                         end do
00625 
00626                         !unk2
00627                         t=2
00628                         do myT=1,endIndex(t)
00629 
00630                                 ! center values
00631                                 UC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00632 
00633                                 !print*,TC,i,j
00634 
00635                                 RHS(t)=UC-0.d0
00636 
00637                         end do
00638 
00639                         !unk3
00640                         t=3
00641                         do myT=1,endIndex(t)
00642 
00643                                 ! center values
00644                                 VC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00645 
00646                                 !print*,TC,i,j
00647 
00648                                 RHS(t)=VC-0.d0
00649 
00650                         end do
00651 
00652 
00653                         do myT=1,maxT
00654                                 
00655                                 t=tVec(myT)
00656 
00657                                 !print*,IM(i,j,t,relative)
00658 
00659                                 discrete(IM(i,j,t,relative))=rhs(t)
00660 
00661 
00662                         end do
00663 
00664                 end do
00665 
00666 
00667         end do
00668         !---------------------------------------------------------------------------------------------
00669         !stop
00670         !---------------------------------------------------------------------------------------------
00671 
00672         !bottom
00673         !bc=4
00674         myBC=4 
00675 
00676         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00677         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00678 
00679 
00680         do myJ=1,maxJ
00681                 
00682                 j=jVec(myJ)
00683                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00684                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00685 
00686                 RHS=0
00687                 do myI=1,maxI
00688 
00689                         i=iVec(myI)
00690                         
00691                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00692                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00693                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00694 
00695                         !unk1
00696                         t=1
00697                         do myT=1,endIndex(t)
00698 
00699                                 ! center values
00700                                 PC=x(IM(i,j,t,absolute))
00701                                 PN=x(IM(i,j+1,t,absolute))
00702                                 PNN=x(IM(i,j+2,t,absolute))
00703                                 
00704                                 VC=x(IM(i,j,t+2,absolute))
00705                                 VN=x(IM(i,j+1,t+2,absolute))
00706                                 VNN=x(IM(i,j+2,t+2,absolute))
00707                                 VNNN=x(IM(i,j+3,t+2,absolute))
00708 
00709                                 !print*,TC,i,j
00710 
00711                                 !RHS(t)=VORC+(UN-UC)/dy
00712 
00713                                 !RHS(t)=VORC+phyOrder*(-3.d0*UC+4.d0*UN-x(IM(i,j+2,t-2,absolute)))/(2.d0*dy)+(1.d0-phyOrder)*(UN-UC)/dy
00714 
00715                                 !RHS(t)=VORC+(-3.d0*UC+4.d0*UN-x(IM(i,j+2,t-2,absolute)))/(2.d0*dy)
00716 
00717                                 !RHS(t)=phyOrder*(Re*(-3.d0*PC+4.d0*PN-1.d0*PNN)/(2.d0*dy)-(VC-5.d0*VN+4.d0*VNN-1.d0*VNNN)/dy**2.d0) +(1.d0-phyOrder)*(PN-PC)/dy ! Re*del[p]/del[y]=del^2v/dely^2
00718                                 !RHS(t)=phyOrder*(Re*(-3.d0*PC+4.d0*PN-1.d0*PNN)/(2.d0*dy)-(VC-5.d0*VN+4.d0*VNN-1.d0*VNNN)/dy**2.d0) +(1.d0-phyOrder)*(Re*(-3.d0*PC+4.d0*PN)/(2.d0*dy)-(VC-5.d0*VN)/dy**2.d0)
00719                                 RHS(t)=phyOrder*(Re*(-3.d0*PC+4.d0*PN-1.d0*PNN)/(2.d0*dy)-(VC-5.d0*VN+4.d0*VNN-1.d0*VNNN)/dy**2.d0) +(1.d0-phyOrder)*(Re*(-3.d0*PC)/(2.d0*dy))
00720                                 RHS(t)=(Re*(-3.d0*PC+4.d0*PN-1.d0*PNN)/(2.d0*dy)-(VC-5.d0*VN+4.d0*VNN-1.d0*VNNN)/dy**2.d0)
00721 
00722 
00723                         end do
00724 
00725                         !unk2
00726                         t=2
00727                         do myT=1,endIndex(t)
00728 
00729                                 ! center values
00730                                 UC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00731 
00732                                 !print*,TC,i,j
00733 
00734                                 RHS(t)=UC-0.d0
00735 
00736                         end do
00737 
00738                         !unk3
00739                         t=3
00740                         do myT=1,endIndex(t)
00741 
00742                                 ! center values
00743                                 VC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00744 
00745                                 !print*,TC,i,j
00746 
00747                                 RHS(t)=VC-0.d0
00748 
00749                         end do
00750 
00751 
00752                         do myT=1,maxT
00753                                 
00754                                 t=tVec(myT)
00755 
00756                                 !print*,IM(i,j,t,relative)
00757 
00758                                 discrete(IM(i,j,t,relative))=rhs(t)
00759 
00760 
00761                         end do
00762 
00763                 end do
00764 
00765 
00766         end do
00767         !---------------------------------------------------------------------------------------------
00768         !stop
00769         !---------------------------------------------------------------------------------------------
00770 
00771         !top
00772         !bc=5
00773         myBC=5 
00774 
00775         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00776         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00777 
00778 
00779         do myJ=1,maxJ
00780                 
00781                 j=jVec(myJ)
00782                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00783                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00784 
00785                 RHS=0
00786                 do myI=1,maxI
00787 
00788                         i=iVec(myI)
00789                         
00790                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00791                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00792                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00793 
00794                         !unk1
00795                         t=1
00796                         do myT=1,endIndex(t)
00797 
00798                                 ! center values
00799                                 PC=x(IM(i,j,t,absolute))
00800                                 PS=x(IM(i,j-1,t,absolute))
00801                                 PSS=x(IM(i,j-1,t,absolute))
00802                                 
00803                                 VC=x(IM(i,j,t+2,absolute))
00804                                 VS=x(IM(i,j-1,t+2,absolute))
00805                                 VSS=x(IM(i,j-2,t+2,absolute))
00806                                 VSSS=x(IM(i,j-3,t+2,absolute))
00807 
00808                                 !print*,TC,i,j
00809 
00810                                 
00811                                 !print*,STFC,STFS,dy,2.d0*(STFC-STFS-dy)/dy**2.d0
00812                                 !stop
00813 
00814                                 !RHS(t)=VORC+(UC-US)/dy
00815                                 !RHS(t)=VORC+phyOrder*(3.d0*UC-4.d0*US+x(IM(i,j-2,t-2,absolute)))/(2.d0*dy)+(1.d0-phyOrder)*(UC-US)/dy
00816                                 !RHS(t)=VORC+(3.d0*UC-4.d0*US+x(IM(i,j-2,t-2,absolute)))/(2.d0*dy)
00817 
00818                                 !RHS(t)=phyOrder*(Re*(3.d0*PC-4.d0*PS+1.d0*PSS)/(2.d0*dy)-(2.d0*VC-5.d0*VS+4.d0*VSS-1.d0*VSSS)/dy**2.d0)+(1.d0-phyOrder)*(PC-PS)/dy ! Re*del[p]/del[y]=del^2v/dely^2
00819                                 !RHS(t)=phyOrder*(Re*(3.d0*PC-4.d0*PS+1.d0*PSS)/(2.d0*dy)-(2.d0*VC-5.d0*VS+4.d0*VSS-1.d0*VSSS)/dy**2.d0)+(1.d0-phyOrder)*(Re*(3.d0*PC-4.d0*PS)/(2.d0*dy)-(2.d0*VC-5.d0*VS)/dy**2.d0)
00820                                 !RHS(t)=phyOrder*(Re*(3.d0*PC-4.d0*PS+1.d0*PSS)/(2.d0*dy)-(2.d0*VC-5.d0*VS+4.d0*VSS-1.d0*VSSS)/dy**2.d0)+(1.d0-phyOrder)*(Re*(3.d0*PC)/(2.d0*dy))
00821                                 RHS(t)=(Re*(3.d0*PC-4.d0*PS+1.d0*PSS)/(2.d0*dy)-(2.d0*VC-5.d0*VS+4.d0*VSS-1.d0*VSSS)/dy**2.d0)
00822 
00823 
00824                         end do
00825 
00826 
00827                         !unk2
00828                         t=2
00829                         do myT=1,endIndex(t)
00830 
00831                                 ! center values
00832                                 UC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00833 
00834                                 !print*,TC,i,j
00835 
00836                                 RHS(t)=UC-1.d0
00837 
00838                         end do
00839 
00840                         !unk3
00841                         t=3
00842                         do myT=1,endIndex(t)
00843 
00844                                 ! center values
00845                                 VC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00846 
00847                                 !print*,TC,i,j
00848 
00849                                 RHS(t)=VC-0.d0
00850 
00851                         end do
00852 
00853 
00854                         do myT=1,maxT
00855                                 
00856                                 t=tVec(myT)
00857 
00858                                 !print*,IM(i,j,t,relative)
00859 
00860                                 discrete(IM(i,j,t,relative))=rhs(t)
00861 
00862 
00863                         end do
00864 
00865                 end do
00866 
00867 
00868         end do
00869         !---------------------------------------------------------------------------------------------
00870         !stop
00871         !---------------------------------------------------------------------------------------------
00872         !bottom corners
00873         !bc=6
00874         myBC=6 
00875 
00876         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00877         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00878 
00879 
00880         do myJ=1,maxJ
00881                 
00882                 j=jVec(myJ)
00883                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00884                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00885 
00886                 RHS=0
00887                 do myI=1,maxI
00888 
00889                         i=iVec(myI)
00890                         
00891                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00892                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00893                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00894 
00895                         !unk1
00896                         t=1
00897                         do myT=1,endIndex(t)
00898 
00899                                 ! center values
00900                                 PC=x(IM(i,j,t,absolute))
00901 
00902                                 RHS(t)=PC
00903 
00904                         end do
00905 
00906 
00907                         !unk1
00908                         t=2
00909                         do myT=1,endIndex(t)
00910 
00911                                 ! center values
00912                                 UC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00913 
00914                                 RHS(t)=UC-0.d0
00915 
00916                         end do
00917 
00918                         !unk1
00919                         t=3
00920                         do myT=1,endIndex(t)
00921 
00922                                 ! center values
00923                                 VC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00924 
00925                                 RHS(t)=VC-0.d0
00926 
00927                         end do
00928 
00929                         do myT=1,maxT
00930                                 
00931                                 t=tVec(myT)
00932 
00933                                 !print*,IM(i,j,t,relative)
00934 
00935                                 discrete(IM(i,j,t,relative))=rhs(t)
00936 
00937 
00938                         end do
00939 
00940                 end do
00941 
00942 
00943         end do
00944         !---------------------------------------------------------------------------------------------
00945         !stop
00946         !---------------------------------------------------------------------------------------------
00947         !top corners
00948         !bc=7
00949         myBC=7 
00950 
00951         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00952         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00953 
00954 
00955         do myJ=1,maxJ
00956                 
00957                 j=jVec(myJ)
00958                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00959                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00960 
00961                 RHS=0
00962                 do myI=1,maxI
00963 
00964                         i=iVec(myI)
00965                         
00966                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00967                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00968                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00969 
00970                         !unk1
00971                         t=1
00972                         do myT=1,endIndex(t)
00973 
00974                                 ! center values
00975                                 PC=x(IM(i,j,t,absolute))
00976 
00977                                 RHS(t)=PC
00978 
00979                         end do
00980 
00981 
00982                         !unk1
00983                         t=2
00984                         do myT=1,endIndex(t)
00985 
00986                                 ! center values
00987                                 UC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00988 
00989                                 RHS(t)=UC-1.d0
00990 
00991                         end do
00992 
00993                         !unk1
00994                         t=3
00995                         do myT=1,endIndex(t)
00996 
00997                                 ! center values
00998                                 VC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00999 
01000                                 RHS(t)=VC-0.d0
01001 
01002                         end do
01003 
01004 
01005                         do myT=1,maxT
01006                                 
01007                                 t=tVec(myT)
01008 
01009                                 !print*,IM(i,j,t,relative)
01010 
01011                                 discrete(IM(i,j,t,relative))=rhs(t)
01012 
01013 
01014                         end do
01015 
01016                 end do
01017 
01018 
01019         end do
01020         !---------------------------------------------------------------------------------------------
01021         !stop
01022 
01023 
01024         end do ! domain 1 loop
01025 
01026 
01027 
01028 
01029 end do ! physics loop
01030 
01031 !#####################################################################
01032 
01033 
01034 deallocate(RHS)
01035 
01036 
01037 
01038 
01039 
01040 
01041 !===================================================================
01042 discreteY=>y
01043 
01044 end function discrete
01045 
01046 
01047 
01048 
01049 
01050 
01051 
01052 
01053 
01054 
01055 
01056 
01057 end module model
D:/demona/2.0/manual/doc/2-31set/problems/ldc_fd/1/model.f90
