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

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) STFc,STFw,STFe,STFs,STFn,stfcpto
00053 real(8) VORc,VORw,VORe,VORs,VORn,vorcpto
00054 
00055 real(8) VORCOREAL
00056 
00057 real(8) dSTFdX,dSTFdY
00058 real(8) dVORdX,dVORdY
00059 real(8) d2STFdX2,d2STFdY2
00060 real(8) d2VORdX2,d2VORdY2
00061 
00062 real(8) dvordtau,dstfdtau
00063 
00064 
00065 real(8) UdVORdX,VdVORdY
00066 
00067 
00068 !real(8) dtau
00069 
00070 real(8) upu,upv
00071 
00072 real(8) myU,myV
00073 
00074 
00075 
00076 real(8) sourceT
00077 real(8) sourceU
00078 
00079 real(8) dx,dy
00080 
00081 real(8),pointer :: myLx,myLy
00082 
00083 real(8),pointer :: Re
00084 
00085 
00086 !< if x+eps*v is to be evaluated, then tempy is updated before calling discrete(1)
00087 !< such that tempy(setVec)=tempy(setVec)+lambda*deltax
00088 !< if discrete(0) is called, original y is evaluated
00089 do i=1,a
00090         discreteY=>tempy
00091 end do
00092 
00093 
00094 
00095 
00096 
00097 
00098 
00099 
00100 discrete=0.d0
00101 
00102 
00103 x=>discreteY
00104 
00105 xPTO=>yPTO
00106 
00107 
00108 discrete=x(setVEc)
00109 
00110 
00111 
00112 ! In this problem, discretization for 1 and 2 is the same although they are located in different
00113 ! part of the domains of different physics
00114 ! they might be coded in a loop so double explicit definition for P1 and P2 is avoided
00115 ! this is not the case for BC 3 since it is defined explicitly
00116 ! however, if and interpolation interface is used, than the bourden is passed to the interface routine
00117 ! so that the discretizations looks the same. After then, P1 and P2 can be defined in a loop.
00118 
00119 
00120 
00121 
00122 
00123 
00124 
00125 
00126 !!!!!!!!!!!!!!
00127 !< assume that  physics-1 starts here
00128 !< also that the domain-1 of physics-1 starts here, as well
00129 !#############################################################################################
00130 !#############################################################################################
00131 !#############################################################################################
00132 !
00133 !############                                   ###
00134 !#############                                  ###
00135 !###            ###                                     ###
00136 !###            ###                                     ###
00137 !###            ###                                     ###
00138 !###            ###                                     ###
00139 !#############          ####            ###
00140 !############           ####            ###
00141 !###                                                    ###
00142 !###                                                    ###
00143 !###                                                    ###
00144 !###                                                    ###
00145 !###                                                    ###
00146 !###                                                    ###     
00147 !
00148 !#############################################################################################
00149 !#############################################################################################
00150 !#############################################################################################
00151 
00152 
00153 !modelPhysics=currentPhysics
00154 !modelDomain=currentDomain
00155 
00156 modelPhysics=1
00157 
00158 
00159 !< the idea is  mymodel%discrete%set(1)
00160 !<                              mymodel%discrete%set(2)
00161 !<                              mymodel%discrete%set(3)
00162 ! with set vec pointing to the set of the physics(..)%domain%(..)%set
00163 
00164 
00165 RE=>myModel%physics(modelPhysics)%physicalParameters(1)
00166 
00167 
00168 
00169 !Re=100.d0
00170 !Re=1000.d0
00171 !Re=2000.d0
00172 !Re=3000.d0
00173 !Re=3500.d0
00174 !Re=5000.d0
00175 !Re=10000.d0
00176 !stop
00177 
00178 
00179 
00180 !print*,'dx,dy=',dx,dy
00181 !stop
00182 
00183 
00184 dof=>myModel%physics(modelPhysics)%dof
00185 
00186 
00187 
00188 
00189 allocate(RHS(dof))
00190 RHS=0.d0
00191 
00192 !#####################################################################
00193 !user supplied variables
00194 !#####################################################################
00195 
00196 do myP=1,myModel%endIndex(modelPhysics)
00197 
00198         !#############################################################################################
00199         !
00200         !########                       ##
00201         !##             ##                      ##
00202         !##             ##                      ##
00203         !##             ##      ###             ##
00204         !##             ##                      ##
00205         !##             ##                      ##
00206         !########                       ##
00207         !
00208         !#############################################################################################
00209         
00210         modelDomain=1
00211 
00212         IM=>myModel%physics(modelPhysics)%domain(modelDomain)%dofMatrix
00213 
00214         myLx=>myModel%physics(modelPhysics)%domain(modelDomain)%Lx
00215         myLy=>myModel%physics(modelPhysics)%domain(modelDomain)%Ly
00216 
00217         dx=myLx/(myModel%physics(modelPhysics)%domain(modelDomain)%maxI-myModel%physics(modelPhysics)%domain(modelDomain)%minI)
00218         dy=myLy/(myModel%physics(modelPhysics)%domain(modelDomain)%maxJ-myModel%physics(modelPhysics)%domain(modelDomain)%minJ)
00219         ! normally, especially in finite volume coordinates should be read so that face lengths and center-to-center distances
00220         ! are calculated. for this uniform finite difference grid dx,dy is enough to proceed
00221 
00222         ! caution, how is dx dy in finite volume!!! was lx for the whole domain with ghostcells 
00223         ! or just the real interior domain w/o ghost cells!!! - if so number of ghostcells connection!!
00224 
00225         !print*,dx,dy
00226 
00227 
00228 
00229         do myD=1,myModel%physics(modelPhysics)%endIndex(modelDomain)
00230 
00231 
00232         !boundary
00233 
00234         !---------------------------------------------------------------------------------------------
00235         !left
00236         !bc=2
00237         myBC=2   
00238 
00239         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00240         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00241 
00242 
00243         do myJ=1,maxJ
00244                 
00245                 j=jVec(myJ)
00246                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00247                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00248 
00249                 RHS=0.d0
00250                 do myI=1,maxI
00251 
00252                         i=iVec(myI)
00253                         
00254                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00255                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00256                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00257 
00258                         !unk1
00259                         t=1
00260                         do myT=1,endIndex(t)
00261 
00262                                 ! center values
00263                                 STFC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00264 
00265                                 !print*,TC,i,j
00266 
00267                                 RHS(t)=STFC-0.d0
00268 
00269                         end do
00270 
00271                         !unk1
00272                         t=2
00273                         do myT=1,endIndex(t)
00274 
00275                                 ! center values
00276                                 STFC=x(IM(i,j,t-1,absolute))
00277                                 VORC=x(IM(i,j,t,absolute))
00278                                 
00279                                 STFE=x(IM(i+1,j,t-1,absolute))
00280 
00281                                 !print*,TC,i,j
00282 
00283                                 RHS(t)=VORC-2.d0*(STFC-STFE)/dx**2.d0
00284 
00285                         end do
00286 
00287 
00288                         do myT=1,maxT
00289                                 
00290                                 t=tVec(myT)
00291 
00292                                 !print*,IM(i,j,t,relative)
00293 
00294                                 discrete(IM(i,j,t,relative))=rhs(t)
00295 
00296                         end do
00297 
00298                 end do
00299 
00300 
00301 
00302         end do
00303         !---------------------------------------------------------------------------------------------
00304         !stop
00305         !---------------------------------------------------------------------------------------------
00306         !right
00307         !bc=3
00308         myBC=3 
00309 
00310         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00311         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00312 
00313 
00314         do myJ=1,maxJ
00315                 
00316                 j=jVec(myJ)
00317                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00318                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00319 
00320                 RHS=0.d0
00321                 do myI=1,maxI
00322 
00323                         i=iVec(myI)
00324                         
00325                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00326                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00327                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00328 
00329                         !unk1
00330                         t=1
00331                         do myT=1,endIndex(t)
00332 
00333                                 ! center values
00334                                 STFC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00335 
00336                                 !print*,TC,i,j
00337 
00338                                 RHS(t)=STFC-0.d0
00339 
00340                         end do
00341 
00342                         !unk1
00343                         t=2
00344                         do myT=1,endIndex(t)
00345 
00346                                 ! center values
00347                                 STFC=x(IM(i,j,t-1,absolute))
00348                                 VORC=x(IM(i,j,t,absolute))
00349                                 
00350                                 STFW=x(IM(i-1,j,t-1,absolute))
00351 
00352                                 !print*,TC,i,j
00353 
00354                                 RHS(t)=VORC-2.d0*(STFC-STFW)/dx**2.d0
00355 
00356                         end do
00357 
00358 
00359                         do myT=1,maxT
00360                                 
00361                                 t=tVec(myT)
00362 
00363                                 !print*,IM(i,j,t,relative)
00364 
00365                                 discrete(IM(i,j,t,relative))=rhs(t)
00366 
00367 
00368                         end do
00369 
00370                 end do
00371 
00372 
00373         end do
00374         !---------------------------------------------------------------------------------------------
00375         !stop
00376         !---------------------------------------------------------------------------------------------
00377 
00378         !bottom
00379         !bc=4
00380         myBC=4 
00381 
00382         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00383         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00384 
00385 
00386         do myJ=1,maxJ
00387                 
00388                 j=jVec(myJ)
00389                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00390                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00391 
00392                 RHS=0.d0
00393                 do myI=1,maxI
00394 
00395                         i=iVec(myI)
00396                         
00397                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00398                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00399                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00400 
00401                         !unk1
00402                         t=1
00403                         do myT=1,endIndex(t)
00404 
00405                                 ! center values
00406                                 STFC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00407 
00408                                 !print*,TC,i,j
00409 
00410                                 RHS(t)=STFC-0.d0
00411 
00412                         end do
00413 
00414                         !unk1
00415                         t=2
00416                         do myT=1,endIndex(t)
00417 
00418                                 ! center values
00419                                 STFC=x(IM(i,j,t-1,absolute))
00420                                 VORC=x(IM(i,j,t,absolute))
00421                                 
00422                                 STFN=x(IM(i,j+1,t-1,absolute))
00423 
00424                                 !print*,TC,i,j
00425 
00426                                 RHS(t)=VORC-2.d0*(STFC-STFN)/dy**2.d0
00427 
00428                         end do
00429 
00430 
00431                         do myT=1,maxT
00432                                 
00433                                 t=tVec(myT)
00434 
00435                                 !print*,IM(i,j,t,relative)
00436 
00437                                 discrete(IM(i,j,t,relative))=rhs(t)
00438 
00439 
00440                         end do
00441 
00442                 end do
00443 
00444 
00445         end do
00446         !---------------------------------------------------------------------------------------------
00447         !stop
00448         !---------------------------------------------------------------------------------------------
00449 
00450         !top
00451         !bc=5
00452         myBC=5 
00453 
00454         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00455         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00456 
00457 
00458         do myJ=1,maxJ
00459                 
00460                 j=jVec(myJ)
00461                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00462                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00463 
00464                 RHS=0.d0
00465                 do myI=1,maxI
00466 
00467                         i=iVec(myI)
00468                         
00469                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00470                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00471                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00472 
00473                         !unk1
00474                         t=1
00475                         do myT=1,endIndex(t)
00476 
00477                                 ! center values
00478                                 STFC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00479 
00480                                 !print*,TC,i,j
00481 
00482                                 RHS(t)=STFC-0.d0
00483 
00484                         end do
00485 
00486                         !unk1
00487                         t=2
00488                         do myT=1,endIndex(t)
00489 
00490                                 ! center values
00491                                 STFC=x(IM(i,j,t-1,absolute))
00492                                 VORC=x(IM(i,j,t,absolute))
00493                                 
00494                                 STFS=x(IM(i,j-1,t-1,absolute))
00495 
00496                                 !print*,TC,i,j
00497 
00498                                 
00499                                 !print*,STFC,STFS,dy,2.d0*(STFC-STFS-dy)/dy**2.d0
00500                                 !stop
00501 
00502                                 RHS(t)=VORC-2.d0*(STFC-STFS-dy)/dy**2.d0
00503 
00504                         end do
00505 
00506 
00507                         do myT=1,maxT
00508                                 
00509                                 t=tVec(myT)
00510 
00511                                 !print*,IM(i,j,t,relative)
00512 
00513                                 discrete(IM(i,j,t,relative))=rhs(t)
00514 
00515 
00516                         end do
00517 
00518                 end do
00519 
00520 
00521         end do
00522         !---------------------------------------------------------------------------------------------
00523         !stop
00524         !---------------------------------------------------------------------------------------------
00525         !corners
00526         !bc=6
00527         myBC=6 
00528 
00529         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00530         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00531 
00532 
00533         do myJ=1,maxJ
00534                 
00535                 j=jVec(myJ)
00536                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00537                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00538 
00539                 RHS=0.d0
00540                 do myI=1,maxI
00541 
00542                         i=iVec(myI)
00543                         
00544                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00545                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00546                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00547 
00548                         !unk1
00549                         t=1
00550                         do myT=1,endIndex(t)
00551 
00552                                 ! center values
00553                                 STFC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00554 
00555                                 RHS(t)=STFC-0.d0
00556 
00557                         end do
00558 
00559                         !unk1
00560                         t=2
00561                         do myT=1,endIndex(t)
00562 
00563                                 ! center values
00564                                 VORC=x(IM(i,j,t,absolute))
00565 
00566                                 RHS(t)=VORC-0.d0
00567 
00568                         end do
00569 
00570 
00571                         do myT=1,maxT
00572                                 
00573                                 t=tVec(myT)
00574 
00575                                 !print*,IM(i,j,t,relative)
00576 
00577                                 discrete(IM(i,j,t,relative))=rhs(t)
00578 
00579 
00580                         end do
00581 
00582                 end do
00583 
00584 
00585         end do
00586         !---------------------------------------------------------------------------------------------
00587         !stop
00588         !---------------------------------------------------------------------------------------------
00589         !interior
00590         !bc=1
00591         myBC=1
00592 
00593         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00594         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00595 
00596         do myJ=1,maxJ
00597                 
00598                 j=jVec(myJ)
00599                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00600                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00601 
00602                 RHS=0.d0
00603                 do myI=1,maxI
00604 
00605                         i=iVec(myI)
00606                         
00607                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00608                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00609                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00610 
00611                         !unk1
00612                         t=1
00613                         do myT=1,endIndex(t)
00614 
00615 
00616                                 STFC=x(IM(i,j,t,absolute))
00617                                 STFW=x(IM(i-1,j,t,absolute))
00618                                 STFE=x(IM(i+1,j,t,absolute))
00619                                 STFS=x(IM(i,j-1,t,absolute))
00620                                 STFN=x(IM(i,j+1,t,absolute))
00621 
00622                                 VORC=x(IM(i,j,t+1,absolute))
00623 
00624                                 stfcpto=xPTO(IM(i,j,t,absolute))
00625 
00626                                 dstfdtau=-0.d0*(stfc-stfcpto)/dtau
00627 
00628                                 d2STFdX2=(STFW+STFE-2.d0*STFC)/(dx**2.d0)
00629 
00630                                 d2STFdY2=(STFN+STFS-2.d0*STFC)/(dy**2.d0)
00631 
00632                                 RHS(t)=dstfdtau+d2STFdX2+d2STFdY2+VORC
00633 
00634 
00635 
00636 
00637 
00638                         end do
00639 
00640                         !unk1
00641                         t=2
00642                         do myT=1,endIndex(t)
00643 
00644                                 STFC=x(IM(i,j,t-1,absolute))
00645                                 STFW=x(IM(i-1,j,t-1,absolute))
00646                                 STFE=x(IM(i+1,j,t-1,absolute))
00647                                 STFS=x(IM(i,j-1,t-1,absolute))
00648                                 STFN=x(IM(i,j+1,t-1,absolute))
00649 
00650 
00651                                 VORC=x(IM(i,j,t,absolute))
00652                                 VORW=x(IM(i-1,j,t,absolute))
00653                                 VORE=x(IM(i+1,j,t,absolute))
00654                                 VORS=x(IM(i,j-1,t,absolute))
00655                                 VORN=x(IM(i,j+1,t,absolute))
00656                                 
00657                                 VORCpto=xpto(IM(i,j,t,absolute))
00658 
00659 
00660                                 dvordtau=0.d0*(vorc-vorcpto)/dtau
00661 
00662                                 dSTFdX=(STFE-STFW)/(2.d0*dx)
00663 
00664                                 dSTFdY=(STFN-STFS)/(2.d0*dy)
00665 
00666                                 !print*,dSTFdX,dSTFdY
00667                                 !stop
00668 
00669                                 !dSTFdX=(STFE-STFC)/(1.d0*dx)
00670 
00671                                 !dSTFdY=(STFC-STFS)/(1.d0*dy)
00672 
00673                                 
00674                                 myU=dSTFdY
00675                                 myV=-dSTFdX
00676 
00677                                 upu=signum(myU)
00678                                 upv=signum(myV)
00679 
00680 
00681                                 UdVORdX=myU*(1.d0-upu)*(VORe-VORc)/(2.d0*dx)+myU*(1.d0+upu)*(VORc-VORw)/(2.d0*dx)
00682                                 VdVORdY=myV*(1.d0-upv)*(VORN-VORc)/(2.d0*dy)+myV*(1.d0+upv)*(VORc-VORS)/(2.d0*dy)
00683 
00684                                 
00685 
00686 
00687                                 dVORdX=(VORE-VORW)/(2.d0*dx)
00688 
00689                                 dVORdY=(VORN-VORS)/(2.d0*dy)
00690 
00691                                 !dVORdX=(VORC-VORW)/(1.d0*dx)
00692 
00693                                 !dVORdY=(VORN-VORC)/(1.d0*dy)
00694 
00695 
00696                                 d2VORdX2=(VORW+VORE-2.d0*VORC)/(dx**2.d0)
00697 
00698                                 d2VORdY2=(VORN+VORS-2.d0*VORC)/(dy**2.d0)
00699 
00700 
00701                                 !RHS(t)=RE*(dvordtau+dSTFdY*dVORdX-dSTFdx*dVORdY)-(d2VORdX2+d2VORdY2)
00702                                 RHS(t)=dvordtau+RE*(dSTFdY*dVORdX-dSTFdx*dVORdY)-(d2VORdX2+d2VORdY2)
00703                                 !RHS(t)=RE*(UdVORdX+VdVORdY)-(d2VORdX2+d2VORdY2)
00704 
00705 
00706                         end do
00707 
00708 
00709 
00710                         do myT=1,maxT
00711 
00712                                 t=tVec(myT)
00713 
00714                                 !print*,IM(i,j,t,relative),rhs(t)
00715 
00716                                 discrete(IM(i,j,t,relative))=rhs(t)
00717 
00718                         end do
00719 
00720                 end do
00721 
00722 
00723         end do
00724         !---------------------------------------------------------------------------------------------
00725         !stop
00726         !---------------------------------------------------------------------------------------------
00727         !artificial left
00728         !bc=7
00729         myBC=7   
00730 
00731         maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
00732         jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec
00733 
00734 
00735         do myJ=1,maxJ
00736                 
00737                 j=jVec(myJ)
00738                 maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
00739                 iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec
00740 
00741                 RHS=0.d0
00742                 do myI=1,maxI
00743 
00744                         i=iVec(myI)
00745                         
00746                         maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
00747                         tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
00748                         endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex
00749 
00750                         !unk1
00751                         t=1
00752                         do myT=1,endIndex(t)
00753 
00754                                 ! center values
00755                                 STFC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))
00756 
00757                                 !print*,TC,i,j
00758 
00759                                 RHS(t)=STFC-0.d0
00760 
00761                         end do
00762 
00763                         !unk1
00764                         t=2
00765                         do myT=1,endIndex(t)
00766 
00767                                 ! center values
00768 
00769                                 VORC=x(IM(i,j,t,absolute))
00770 
00771                                 !print*,TC,i,j
00772 
00773                                 RHS(t)=VORC-0.d0
00774 
00775                         end do
00776 
00777 
00778                         do myT=1,maxT
00779                                 
00780                                 t=tVec(myT)
00781 
00782                                 !print*,IM(i,j,t,relative)
00783 
00784                                 discrete(IM(i,j,t,relative))=rhs(t)
00785 
00786                         end do
00787 
00788                 end do
00789 
00790 
00791 
00792         end do
00793         !---------------------------------------------------------------------------------------------
00794         !stop
00795 
00796         end do ! domain 1 loop
00797 
00798 
00799 
00800 
00801 end do ! physics loop
00802 
00803 !#####################################################################
00804 
00805 
00806 deallocate(RHS)
00807 
00808 
00809 
00810 
00811 
00812 
00813 !===================================================================
00814 discreteY=>y
00815 
00816 end function discrete
00817 
00818 
00819 
00820 
00821 
00822 
00823 
00824 
00825 
00826 
00827 
00828 
00829 end module model
D:/demona/2.0/manual/doc/2-31set/problems/sf_vor/1/model.f90
