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

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