!======================================================================== ! 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 problem parameters ! qinit - Initialize field variables ! physflux - Compute the physical flux ! ! Comments: This module is problem specific and must be written by the ! user. For examples 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 equation, constant velocity advection of an ! initial perturbation ! 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(3) :: 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) READ(7,*) ubar(3) 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,k REAL (KIND=xPrec) :: x,y,z ! Piecewise constant with ! q = 1.0 if 12 < x < 20 and 12 < y < 20 and 12 < z < 20 ! 0.1 otherwise z=Info%Xlower(3)+Info%dX(3)/2.0 DO k=1,Info%mX(3) 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>12 .AND. x<20 .AND. y>12 .AND. y<20 .AND. z>12 .AND. z<20) THEN Info%q(i,j,k,1,1)=1. ELSE Info%q(i,j,k,1,1)=0. END IF x=x+Info%dX(1) END DO y=y+Info%dX(2) END DO z=z+Info%dX(3) END DO END SUBROUTINE qinit SUBROUTINE b4step(Info) ! Interface declarations TYPE (NodeInfo) :: Info ! Internal declarations 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 ! Internal declarations END SUBROUTINE setaux SUBROUTINE src(Info) ! Interface declarations TYPE (NodeInfo) :: Info ! Internal declarations END SUBROUTINE src ! ----------------------------------------------------------------------- ! ! Compute physical fluxes, speeds for the scalar equation ! ! q_t + u*q_x + v*q_y + w*q_z = 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 NodeInfoDef.f90 for definitions TYPE (NodeInfo) :: Info ! Data associated with this grid REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: q,f,s REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: A ! Internal declarations INTEGER mbc,mx,nq,n INTEGER, POINTER, DIMENSION(:) :: I ! Current number of interior cells, boundary cells mx=Info%mXnow; mbc=Info%mbc; nq=Info%NrVars ! Index range I=>Info%I1D SELECT CASE (irequest) CASE (Initialize) CASE (Finalize) CASE (RequestFluxes) f(I,1:nq)=ubar(ixy)*q(I,1:nq) CASE (RequestSpeeds) s(I,1:nq)=ubar(ixy) CASE (RequestJacobians) A=0. DO n=1,nq A(I,n,n)=ubar(ixy) END DO END SELECT END SUBROUTINE physflux !***************** END MODULE problem !*****************