!======================================================================== ! BEARCLAW Boundary Embedded Adaptive Refinement Conservation LAW package !======================================================================== ! (c) Copyright Sorin Mitran, 2000 ! Department of Applied Mathematics ! University of Washington ! mitran@amath.washington.edu ! ----------------------------------------------------------------------- ! This code may be freely used for educational and research purposes. ! For any other use please contact the author. ! ----------------------------------------------------------------------- ! File: problem.f90 ! Purpose: Definition of problem specific data and procedures ! Contains: ! ! setprob - Initialize overall problem parameters ! qinit - Initialize field variables on root grids ! b4step - Called prior to each time step ! setaux - Initialize auxilliary variable fields ! src - Compute source terms ! physflux - Compute the physical flux, related quantities ! Compute wave decomposition ! ! Comments: This module is problem specific and must be written by the ! user. For exampels see the BEARCLAW documentation and also ! the examples and applications directories ! ! The encapsulation of all problem specific data and routines ! and data into a module is useful in exposing data that needs ! to be shared. ! Application: Advection of a square ! Contains: ! Revision History: Ver. 1.0 Oct. 2000 Sorin Mitran ! ----------------------------------------------------------------------- !************* MODULE problem !************* ! User definitions of data structures associated with each node USE NodeInfoDef IMPLICIT NONE ! It's safer to require explicit declarations SAVE ! Save module information PRIVATE ! Everything is implicitly private, i.e. accessible only ! to this module. ! These are the public entry points (acessible to other modules) PUBLIC setprob,qinit,b4step,setaux,src,physflux,problemBC,afterstep,problemIO ! Problem specific global variables: ! - Problem parameters REAL (KIND=qPrec), DIMENSION(2) :: ubar !======= CONTAINS !======= SUBROUTINE setprob ! User definitions of data structures associated with each node USE NodeInfoDef IMPLICIT NONE ! Read the velocity for scalar advection problem OPEN(unit=7,file='setprob.data',status='old',form='formatted') READ(7,*) ubar(1) READ(7,*) ubar(2) CLOSE(7) END SUBROUTINE setprob SUBROUTINE problemBC(Info) ! User definitions of data structures associated with each node USE NodeInfoDef ! Interface declarations TYPE (NodeInfo) :: Info ! Data associated with this grid ! Internal declarations END SUBROUTINE problemBC SUBROUTINE qinit(Info) ! User definitions of data structures associated with each node USE NodeInfoDef ! Interface declarations TYPE (NodeInfo) :: Info ! Data associated with this grid ! Internal declarations INTEGER i,j REAL (KIND=xPrec) :: x,y ! ,r ! Piecewise constant with ! q = 1.0 if 12 < x < 20 and 12 < y < 20 ! 0.0 otherwise y=Info%Xlower(2)+Info%dX(2)/2.0 DO j=1,Info%mX(2) x=Info%Xlower(1)+Info%dX(1)/2.0 DO i=1,Info%mX(1) IF ( x>8 .AND. x<24 .AND. y>8 .AND. y<24 ) THEN Info%q(i,j,1,1,1)=1.0 ELSE Info%q(i,j,1,1,1)=0.0 END IF !r=SQRT((x-12)**2+(y-12)**2) !Info%q(i,j,1,1,1)=exp(-1.5*r) x=x+Info%dX(1) END DO y=y+Info%dX(2) END DO END SUBROUTINE qinit SUBROUTINE b4step(Info) ! Interface declarations TYPE (NodeInfo) :: Info END SUBROUTINE b4step SUBROUTINE afterstep(Info) ! Interface declarations TYPE (NodeInfo) :: Info ! Include operations that need to be carried out after each time step here: END SUBROUTINE afterstep SUBROUTINE problemIO(nframe,tnow,IOrequest,Info,qmax,qmin) ! Interface declarations INTEGER :: nframe,IOrequest; REAL (KIND=qPrec) :: tnow TYPE (NodeInfo), OPTIONAL :: Info REAL (KIND=qPrec), OPTIONAL, DIMENSION(:) :: qmax,qmin INTEGER, PARAMETER :: UserBeforeGridIO=-1, UserAfterGridIO=-2 INTEGER, PARAMETER :: UserIOMinMax=0 ! This routine provides for output of computation snapshots in user-defined ! formats. ! IOrequest is a function parameter with values: ! -1 UserBeforeGridIO - indicates routine is being called prior to parsing ! of tree of grids. This is where output files ! should be opened ! -2 UserAfterGridIO - indicates routine is being called after parsing ! of tree of grids. This is where output files ! should be closed ! -3 UserOutputq - indicates routine should output the field ! variables in user specified format ! 0 UserIOMinMax - routine is being called after qmin and qmax have ! been determined. These are the min/max field ! variable values. If output of quantities different ! from the field variables is desired, qmin & ! qmax should be defined with respect to those quantities ! 1,2 ... - routine is being called to output field variable ! 1,2,... on the current grid in conjunction with ! HDF support routines in BEARIO. Only valid for ! UserHDF output style END SUBROUTINE problemIO SUBROUTINE setaux(Info) ! Interface declarations TYPE (NodeInfo) :: Info END SUBROUTINE setaux SUBROUTINE src(Info) ! Interface declarations TYPE (NodeInfo) :: Info END SUBROUTINE src ! ----------------------------------------------------------------------- ! ! Compute the Riemann problem solution for the advection equation ! ! q_t + u*q_x + v*q_y = 0 ! ! ----------------------------------------------------------------------- SUBROUTINE physflux(ixy,indx,irequest,Info,q,f,s,A) ! User definitions of data structures associated with each node USE NodeInfoDef ! Interface declarations INTEGER ixy ! The dimension along which to solve the normal Riemann problem INTEGER indx(MaxDims) ! Indices of this slice, indx(ixy) has no significance, other ! indices give position of this 1D slice in index space INTEGER irequest ! Request code, see NodeInfoGlobal.f90 for definitions ! Additional request codes may be defined in solver modules TYPE (NodeInfo) :: Info ! Data associated with this grid REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: q,f,s REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: A ! Shorter local names for the wave propagation quantities REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: Apdq,Amdq,speed REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: Asdq,BmAsdq,BpAsdq REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: wave ! Internal declarations INTEGER mbc,mx,nq,tdir,j INTEGER, POINTER, DIMENSION(:) :: I REAL (KIND=qPrec) :: stran,stranm,stranp,zero=0. ! Current number of interior cells, boundary cells mx=Info%mXnow; mbc=Info%mbc; nq=Info%NrVars; tdir=Info%tdir ! Index range. Contains indices of interfaces where the ! solution to the Riemann problem is requested between ! state I-1 to the left and I to the right I=>Info%I1D ! Associate local names with Info components Apdq=>Info%Apdq; Amdq=>Info%Amdq; speed=>Info%speed; wave=>Info%wave Asdq=>Info%Asdq; BmAsdq=>Info%BmAsdq; BpAsdq=>Info%BpAsdq SELECT CASE (irequest) CASE (Initialize) CASE (Finalize) CASE (RequestFluxes) f(I,1:nq)=ubar(ixy)*q(I,1:nq) CASE (RequestNormalWaves) Apdq=0.; Amdq=0.; speed=0.; wave=0. wave(I,1,1)=q(I,1)-q(I-1,1) speed(I,1)=ubar(ixy) ! Flux differences WHERE (speed(I,1)<0.) Amdq(I,1)=speed(I,1)*wave(I,1,1) ELSEWHERE Apdq(I,1)=speed(I,1)*wave(I,1,1) END WHERE ! Elegant way to do it, but the NAG f95 v4.0 compiler has a bug ! and leaks memory on the WHERE line ! Leak was caught by memprof, http://people.redhat.com/otaylor/memprof/ ! NAG f95 4.2 does not leak memory. If you suspect a memory leak ! remove comments from code below: !DO j=2-mbc,mx+mbc ! IF (speed(j,1) < 0.) THEN ! Amdq(j,1)=speed(j,1)*wave(j,1,1) ! ELSE ! Apdq(j,1)=speed(j,1)*wave(j,1,1) ! END IF !END DO CASE (RequestTransverseWaves) BmAsdq=0.; BpAsdq=0. stran=ubar(tdir) stranm = MIN(stran, zero) stranp = MAX(stran, zero) BmAsdq(I,1) = stranm * Asdq(I,1) BpAsdq(I,1) = stranp * Asdq(I,1) CASE DEFAULT PRINT *,'Error: the solver module Scheme.f90 emitted a request that' PRINT *,'is not implemented in the problem module problem.f90' STOP END SELECT END SUBROUTINE physflux !***************** END MODULE problem !*****************