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