!> @dir PulsedJet !! @brief Contains files necessary for module PulsedJet !> @file problem.f90 !! @brief Main file for module PulsedJet !> @defgroup PulsedJet PulsedJet Module !! @brief Module for simulating a pulsed 2.5D Jet !! @ingroup Modules !> Pulsed Jet problem module !! @ingroup 1DJet MODULE Problem USE GlobalDeclarations USE PhysicsDeclarations USE DataDeclarations USE CoolingSrc USE ZCooling USE Chemistry USE Refinements USE Shapes USE Fields USE Emissions USE Projections USE ProcessingDeclarations USE Ambients IMPLICIT NONE PUBLIC ProblemModuleInit, ProblemGridInit, & ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep REAL(KIND=qPREC) :: namb, njet, tamb, vjet, Rjet, Bm, vrange, newvjet, period, pamb, alpha, & Betam, Rm, freq, amp, yupper, theta, ycoeff, xzcoeff, precfreq, precperiod, CoolingRes REAL(KIND=qPREC), ALLOCATABLE :: randoms(:), vtimes(:) INTEGER :: Bgeom, nrandoms, seed, vtype, simtype, N INTEGER, PARAMETER :: toroidal = 1, helical = 2 INTEGER, PARAMETER :: constant = 0, sinusoidal = 1, random = 2 INTEGER, PARAMETER :: single = 0, colliding_same = 1, colliding_opp = 2 LOGICAL :: ljet CONTAINS !> Initializes module variables SUBROUTINE ProblemModuleInit() TYPE(AmbientDef), POINTER :: Ambient TYPE(RefinementDef), POINTER :: InjectRefine, TopInjectRefine, CoolingLengthRefinement TYPE(ProjectionDef), POINTER :: Projection REAL(KIND=qPREC) :: rand INTEGER :: i, jj NAMELIST/problemdata/ namb, tamb, njet, vjet, vrange, Rjet, Rm, Betam, Bgeom, vtype, freq, & theta, precfreq, seed, simtype, N, CoolingRes OPEN(PROBLEM_DATA_HANDLE,FILE='problem.data',STATUS='old') READ(PROBLEM_DATA_HANDLE,NML=problemdata) CLOSE(PROBLEM_DATA_HANDLE) ! CALL AddDiagnosticVar(OI_Field, 'OI') ! 6300 ! CALL AddDiagnosticVar(NII_Field, 'NII') ! 6583 ! CALL AddDiagnosticVar(SII_6716_Field, 'SII_6716') ! CALL AddDiagnosticVar(SII_6731_Field, 'SII_6731') ! CALL AddDiagnosticVar(Halpha_Field, 'Halpha') ! CALL AddDiagnosticVar(CellSize_Field, 'dx') ! CALL AddDiagnosticVar(CoolingLength_Field, 'CoolingLength') ! Set up some parameters for the velocity pulses SELECT CASE(vtype) CASE(constant) CASE(sinusoidal) period = final_time/freq amp = vrange/vjet CASE(random) ! set up an array of pseudo-random numbers to use for random jet velocities nrandoms = NINT(freq) period = final_time/REAL(nrandoms,qPrec) ALLOCATE(randoms(nrandoms)) ALLOCATE(vtimes(nrandoms)) CALL RANDOM_SEED(seed) DO jj=1, nrandoms CALL RANDOM_NUMBER(rand) randoms(jj) = rand vtimes(jj) = jj*period END DO END SELECT pamb = namb*Boltzmann*tamb IF(lMHD) THEN Bm = SQRT(8d0*Pi*pamb/Betam)/Bscale alpha = 1d0 - (Rm/Rjet)**2d0/Betam ELSE Bm = 0d0 alpha = 1d0 Rm = -10d0 END IF ! Set up parameters related to a precessing injection velocity theta = theta*Pi/180d0 ycoeff = SQRT(1d0/(1d0 + 2d0*TAN(theta)**2d0)) xzcoeff = SQRT(1d0 - ycoeff**2d0) precperiod = final_time/precfreq ! Set up ambient CALL CreateAmbient(Ambient) Ambient%density= namb/nScale Ambient%pressure= pamb/pScale ! Set up refinement CALL ClearAllRefinements() CALL CreateRefinement(CoolingLengthRefinement) CoolingLengthRefinement%field=CoolingTime_Field CoolingLengthRefinement%limit = LESSTHAN CoolingLengthRefinement%Threshold(0:MaxLevel)=CoolingRes*levels(0:MaxLevel)%dx/(0.1*vjet/VelScale) CALL CreateRefinement(InjectRefine) CALL CreateShape(InjectRefine%Shape) InjectRefine%Shape%type = RECTANGULAR_PRISM InjectRefine%Shape%position = (/ 0d0, 0d0, 0d0 /) InjectRefine%Shape%size_param = (/Rjet, 1*levels(0)%dx, Rjet/) InjectRefine%Maxlevel = 3 CALL SetShapeBounds(InjectRefine%Shape) IF (simtype /= single) THEN yupper = GxBounds(2,2) CALL CreateRefinement(TopInjectRefine) CALL CreateShape(TopInjectRefine%Shape) TopInjectRefine%Shape%type = RECTANGULAR_PRISM TopInjectRefine%Shape%position = (/ 0d0, yupper, 0d0 /) TopInjectRefine%Shape%size_param = Rjet InjectRefine%Maxlevel = 3 CALL SetShapeBounds(TopInjectRefine%Shape) END IF ! Set up projections to make emission maps CALL CreateProjection(Projection) IF (iCylindrical == 0) THEN Projection%dim = 3 ELSE Projection%dim = CYLINDRICAL_PROJECTION END IF Projection%PlotLevel = MIN(MaxLevel, 4) Projection%field(1)%id = Halpha_Field Projection%field(1)%name = "Halpha" Projection%field(2)%id = SII_6716_Field Projection%field(2)%name = "SII_6716" Projection%field(3)%id = SII_6731_Field Projection%field(3)%name = "SII_6731" ! Projection%field(4)%id = NII_Field ! Projection%field(4)%name = "NII" END SUBROUTINE ProblemModuleInit !> Applies initial conditions !! @param Info Info object SUBROUTINE ProblemGridInit(Info) TYPE(InfoDef) :: Info CALL PlaceJet(Info, levels(Info%level)%dx) END SUBROUTINE ProblemGridInit !> Applies Boundary conditions !! @param Info Info object SUBROUTINE ProblemBeforeStep(Info) TYPE(InfoDef) :: Info CALL PlaceJet(Info, 0d0) END SUBROUTINE ProblemBeforeStep SUBROUTINE PlaceJet(Info, yjet) TYPE(InfoDef) :: Info REAL(KIND=qPREC), POINTER, DIMENSION(:,:,:,:) :: q INTEGER :: nn, i, j, k, rmbc, zrmbc, level, mx, my, mz, ii, jj, kk, Nz REAL(KIND=qPREC) :: x, y , z, r, phi, dx, xlower, ylower, zlower, dt, time,& Bx, By, Bz, randv, pjet, Bjet, nH, nHe, Ay, ptot, padd, ddx, xx, yy, zz, rr, yjet, hdx, pB, sumpB, ddz REAL(KIND=qPREC),ALLOCATABLE :: A(:,:,:,:), BzAux(:,:) REAL(KIND=qPREC), DIMENSION(1:NrHydroVars) :: newq, sumq level=Info%level q=>Info%q rmbc = levels(level)%gmbc(levels(level)%step) zrmbc = 0 mx = Info%mX(1) my = Info%mX(2) mz = Info%mX(3) dx = levels(level)%dX ddx = dx/N ddz = 0 Nz = 1 IF (nDim == 3) THEN zrmbc = rmbc Nz=N ddz=ddx END IF xlower = Info%xBounds(1,1) ylower = Info%xBounds(2,1) zlower = Info%xBounds(3,1) hdx=half*dx dt = levels(level)%dt time = levels(level)%tnow ! IF (time <= dt) RETURN ! don't want to overwrite initial conditions phi = 0.0 ! sets new jet velocity depending on time and type of velocity changes SELECT CASE(vtype) CASE(constant) newvjet = vjet CASE(sinusoidal) newvjet = vjet*(1d0 + amp*SIN(2d0*Pi*time/period)) phi = 2d0*Pi*time/precperiod CASE(random) DO nn=1, nrandoms IF (time == vtimes(nn)) THEN randv = ((randoms(nn)*2d0)-1d0)*vrange newvjet = vjet + randv ELSE IF (time < vtimes(1)) THEN newvjet = vjet END IF END DO END SELECT DO j=1-rmbc, my+rmbc y = ylower + (REAL(j,xPrec) - 1d0)*dx IF (simtype /= single) THEN IF (y >= 0d0+yjet .AND. y <= yupper-yjet-dx) cycle ELSE IF (y >= 0d0+yjet) THEN ! write(*,*) j, y cycle END IF END IF DO i=1-rmbc, mx+rmbc x = xlower + (REAL(i,xPrec) - 1d0)*dx DO k=1-zrmbc, mz+zrmbc z = zlower + (REAL(k,xPrec) - 1d0)*dx sumq=0d0 sumpB=0d0 CALL Cons_To_Prim(q(i,j,k,:)) DO jj=1, N yy = y + (REAL(jj,xPrec) - half)*ddx ! ONLY DO FOLLOWING INSIDE THE -y BOUNDARY. ALSO, ! ABOVE THE +y BOUNDARY IF THIS IS A COLLIDING JETS RUN IF (simtype /= single) THEN ljet = (yy <= 0d0+yjet .AND. yy >= yupper-yjet) ELSE ljet = (yy <= 0d0+yjet) ! write(*,'(2I6,2E15.4,L2)') j,jj,y,yy,ljet END IF DO ii=1, N xx = x + (REAL(ii,xPrec) - half)*ddx DO kk=1, Nz IF (ljet) THEN IF (nDim == 2) THEN rr = xx ELSE IF (nDim == 3) THEN zz = z + (REAL(kk,xPrec) - half)*ddz rr = SQRT(xx**2d0+zz**2d0) END IF IF (rr < Rjet) THEN newq(1) = njet/nScale newq(ivx) = xzcoeff*COS(phi)*newvjet/velScale newq(ivy) = ycoeff*newvjet/velScale ! Second jet injection region, if this is a colliding jets run IF (ivz /= 0) THEN IF (simtype /= single .AND. y > yupper-yjet) newq(ivy) = -newq(ivy) newq(ivz) = xzcoeff*SIN(phi)*newvjet/velScale END IF IF (rr < Rm) THEN pjet = pamb*(alpha + 2d0/Betam*(1d0 - (rr/Rm)**2d0))/pScale ELSE pjet = alpha*pamb/pScale END IF newq(iE) = pjet !/gamma1 + half*SUM(newq(ivx:ivz)**2d0)/newq(1) IF (rr < Rm) THEN pB= .5*(Bm*rr/Rm)**2d0 ELSE pB = .5*(Bm*Rm/rr)**2d0 END IF IF (lTrackHydrogen) THEN nH = newq(1)/(1d0+xH+HeFrac*(1d0+xHe+2d0*xHeII)+ZFrac*(1d0+xZ)) newq(iH) = (1d0-xH)*nH newq(iHII) = xH*nH IF (lTrackHelium) THEN nHe = HeFrac*nH newq(iHe) = (1d0-xHe-xHeII)*nHe newq(iHeII) = xHe*nHe newq(iHeIII) = xHeII*nHe END IF END IF ELSE newq=Info%q(i,j,k,1:NrHydroVars) ! newq(1) = newq(1) + Info%q(i,j,k,1,1) !namb/nScale ! newq(ivx:ivz) = newq(ivx:ivz) + 0d0 ! newq(iE) = newq(iE) + pamb/pScale/gamma1 pB=0d0 END IF ELSE newq=Info%q(i,j,k,1:NrHydroVars) pB=0d0 END IF sumq=sumq+newq sumpB=sumpB+pB END DO END DO END DO q(i,j,k,1:NrHydroVars) = sumq/N**nDim CALL Prim_To_Cons(q(i,j,k,:)) q(i,j,k,iE)=q(i,j,k,iE)+sumpB/N**nDim END DO END DO END DO IF (lMHD) THEN ALLOCATE(A(1-rmbc:mx+rmbc+1,1-rmbc:my+rmbc+1,1-zrmbc:mz+zrmbc+1,3)) IF (nDim == 2) ALLOCATE(BzAux(1-rmbc:mx+rmbc+1,1-rmbc:my+rmbc+1)) A(:,:,:,1) = 0d0 A(:,:,:,3) = 0d0 DO j=1-rmbc, my+rmbc+1 y = ylower + (REAL(j,xPrec) - half)*dx ljet = .false. ! Check if inside jet injection region(s) IF (simtype /= single) THEN IF (y < 0d0+yjet .OR. y > yupper-yjet) ljet = .true. ELSE IF (y < 0d0+yjet) ljet = .true. END IF DO i=1-rmbc, mx+rmbc+1 x = xlower + (REAL(i,xPrec) - 1d0)*dx DO k=1-zrmbc, mz+zrmbc+1 ! A M B I E N T Ay = -Bm*Rm*(LOG(Rjet/Rm)+half) IF (ljet) THEN IF (nDim == 2) THEN r = x ELSE IF (nDim == 3) THEN z = zlower + (REAL(k,xPrec) - 1d0)*dx r = SQRT(x**2d0 + z**2d0) END IF IF (r < Rjet) THEN IF (r < Rm) THEN Ay = -half*Bm*r**2d0/Rm ELSE Ay = -Bm*Rm*(LOG(r/Rm)+half) END IF END IF END IF IF (simtype == colliding_opp .AND. y > yupper-dx) Ay = -Ay A(i,j,k,2) = Ay END DO END DO END DO !Then take the curl of A to get B IF (nDim == 2) THEN Info%aux(:,:,:,:) = 0d0 ELSE IF (nDim == 3) THEN DO j=1-rmbc, my+rmbc y = ylower + (REAL(j,xPrec) - half)*dx IF (simtype /= single) THEN IF (y >= 0d0+yjet .AND. y <= yupper-yjet) CYCLE ELSE IF (y >= 0d0+yjet) CYCLE END IF DO i=1-rmbc, mx+rmbc+1 DO k=1-zrmbc, mz+zrmbc Info%aux(i,j,k,1) = (A(i,j+1,k,3)-A(i,j,k,3))/dx - (A(i,j,k+1,2)-A(i,j,k,2))/dx END DO END DO END DO DO j=1-rmbc, my+rmbc+1 y = ylower + (REAL(j,xPrec) - 1d0)*dx IF (simtype /= single) THEN IF (y >= 0d0+yjet .AND. y <= yupper+yjet) CYCLE ELSE IF (y >= 0d0+yjet) CYCLE END IF DO i=1-rmbc, mx+rmbc DO k=1-zrmbc, mz+zrmbc Info%aux(i,j,k,2) = (A(i,j,k+1,1)-A(i,j,k,1))/dx - (A(i+1,j,k,3)-A(i,j,k,3))/dx END DO END DO END DO END IF DO j=1-rmbc, my+rmbc y = ylower + (REAL(j,xPrec) - half)*dx IF (simtype /= single) THEN IF (y >= 0d0+yjet .AND. y <= yupper-yjet) CYCLE ELSE IF (y >= 0d0+yjet) CYCLE END IF DO i=1-rmbc, mx+rmbc DO k=1-zrmbc, mz+zrmbc+1 IF (nDim == 2) THEN BzAux(i,j) = (A(i+1,j,1,2)-A(i,j,1,2))/dx - (A(i,j+1,1,1)-A(i,j,1,1))/dx ELSE IF (nDim == 3) THEN Info%aux(i,j,k,3) = (A(i+1,j,k,2)-A(i,j,k,2))/dx - (A(i,j+1,k,1)-A(i,j,k,1))/dx END IF END DO END DO END DO DEALLOCATE(A) ! Average aux fields to get B-fields for q, update energy DO j=1-rmbc, my+rmbc y = ylower + (REAL(j,xPrec) - half)*dx IF (simtype /= single) THEN IF (y >= 0d0+yjet .AND. y <= yupper-yjet) CYCLE ELSE IF (y >= 0d0+yjet) CYCLE END IF DO i=1-rmbc, mx+rmbc DO k=1-zrmbc, mz+zrmbc IF (nDim == 2) THEN q(i,j,k,iBx:iBy) = 0d0 q(i,j,k,iBz) = BzAux(i,j) ELSE IF (nDim == 3) THEN q(i,j,k,iBx) = half*(Info%aux(i,j,k,1)+Info%aux(i+1,j,k,1)) q(i,j,k,iBy) = half*(Info%aux(i,j,k,2)+Info%aux(i,j+1,k,2)) q(i,j,k,iBz) = half*(Info%aux(i,j,k,3)+Info%aux(i,j,k+1,3)) END IF ! q(i,j,k,iE) = q(i,j,k,iE) + half*SUM(q(i,j,k,iBx:iBz)**2d0) ! Pressure Check ! ptot = (gamma1*(q(i,j,k,iE) - half*SUM(q(i,j,k,2:m_high)**2d0)/q(i,j,k,1)) + & ! gamma13*half*SUM(q(i,j,k,iBx:iBz)**2d0))*pScale ! IF (ptot < pamb) THEN ! padd = pamb - ptot ! q(i,j,k,iE) = q(i,j,k,iE) + padd/gamma1/pScale ! END IF END DO END DO END DO IF (nDim == 2) DEALLOCATE(BzAux) END IF END SUBROUTINE PlaceJet !> Could be used to update grids pre-output !! @param Info Info Object SUBROUTINE ProblemAfterStep(Info) TYPE(InfoDef) :: Info END SUBROUTINE ProblemAfterStep !> Could be used to set force refinement !! @param Info Info object SUBROUTINE ProblemSetErrFlag(Info) TYPE(InfoDef) :: Info END SUBROUTINE ProblemSetErrFlag SUBROUTINE ProblemBeforeGlobalStep(n) INTEGER :: n END SUBROUTINE ProblemBeforeGlobalStep END MODULE Problem