cart-densityfield_fluct.F90

module problem

  ! Module for Capreole2D (F90)
  ! Author: Garrelt Mellema
  ! Date: 2015-08-11 (2004-05-11)

  ! This is the problem module. It contains the routines which define the
  ! problem being solved:

  ! Contents:
  ! init_problem - sets up the hydro variables according to the specified
  !                  problem
  ! inflow - sets the inflow boundary conditions

  ! This version: specify an average density (n) and a field of 
  !    fluctuations (delta). The density field used will be n*(1+delta).
  !    To avoid negative densities we impose a lower density limit eta
  !    so that the minimum density is eta*n.
  !    Temperature can be set isobarically (specify average temperature)
  !    or isothermal. 
  !    Coordinates: cartesian coordinates.

  use file_admin, only: stdinput
  use precision, only: dp
  use cgsconstants, only: m_p
  use my_mpi
  use sizes, only: mbc,neq,RHO,RHVX,RHVY,RHVZ,EN,nrofDim
  use scaling, only: SCDENS, SCMOME, SCENER
  use atomic, only: gamma1
  use abundances, only: mu
  use mesh, only: sx,ex,sy,ey,sz,ez,meshx,meshy,meshz
  !use grid  
  use hydro, only: state,pressr,set_state_pointer,NEW,OLD,restart_state
  use boundary, only: boundaries,REFLECTIVE,OUTFLOW,PROBLEM_DEF,X_IN,X_OUT,Y_IN, &
       Y_OUT,Z_IN,Z_OUT
  use ionic, only: init_ionic
  use tped

  implicit none

  !> lowest density is minimum_fluctuation*average_density
  real(kind=dp),parameter :: minimum_fluctuation=0.1 
  integer, dimension(nrofDim,2) :: domainboundaryconditions

contains

  subroutine init_problem (restart,restartfile)
    
    ! This routine initializes all hydro variables
    
    ! This may be a fresh start or a restart of a saved run
    
    logical,intent(in) :: restart ! tells you whether it's a new run or a restart
    character(len=19),intent(in) :: restartfile

    ! Local variables
    real(kind=dp) :: r_interface ! dummy needed for calling init_ionic
    integer :: ierror

    ! Set domain boundary conditions
    domainboundaryconditions(2:3,:)=REFLECTIVE    
    domainboundaryconditions(1,2)=OUTFLOW
    domainboundaryconditions(1,1)=OUTFLOW

    if (.not.restart) then ! Fresh start

       call fresh_start_state( )
       
    else

       call restart_state(restartfile,ierror)
       state(:,:,:,RHO)=state(:,:,:,RHO)/scdens
       state(:,:,:,RHVX)=state(:,:,:,RHVX)/scmome
       state(:,:,:,RHVY)=state(:,:,:,RHVY)/scmome
       state(:,:,:,RHVZ)=state(:,:,:,RHVZ)/scmome
       state(:,:,:,EN)=state(:,:,:,EN)/scener
       call boundaries(OLD,domainboundaryconditions,problemboundary) ! Fill boundary conditions

    endif
       
    ! Initialize the ionic concentrations
    call init_ionic(restart,r_interface)

  end subroutine init_problem

  !==========================================================================

  subroutine fresh_start_state ( )
    
    ! This routine initializes all hydro variables for a fresh start
    
    ! Case: density fields
    
    ! smoothing parameter for making soft-edged clump
    real(kind=dp),parameter :: eta=0.1_dp
    
    ! e variables are environment
    real(kind=dp)    :: average_temperature
    real(kind=dp)    :: epressr

    integer :: mx1,my1,mz1,temperature_setup
    character(len=512) :: densityfluctuations_file
    real,allocatable,dimension(:,:,:) :: tempdens
    real(kind=dp) :: average_density_cm3
    real(kind=dp) :: maxdens
    integer :: i,j,k,nitt,ieq

#ifdef MPI       
    integer :: status(MPI_STATUS_SIZE)
    integer :: request,nextproc
    integer,parameter :: outputcircle=601
    integer :: ierror
#endif

    ! Ask for the input if you are processor 0.
    if (rank == 0) then
       write (*,'(A,$)') '1) Average density (cm^-3): '
       read (stdinput,*) average_density_cm3
       write (*,'(A,$)') '1) File with density fluctuations: '
       read (stdinput,*) densityfluctuations_file
       write (*,'(A,$)') '2) Temperature setup: '
       write (*,'(A,$)') '   1: isobaric: '
       write (*,'(A,$)') '   2: isothermal: '
       read (stdinput,*) temperature_setup
       if (temperature_setup == 1) then
          write (*,'(A,$)') '3) Average temperature: '
          read (stdinput,*) average_temperature
       else
          write (*,'(A,$)') '3) Temperature: '
          read (stdinput,*) average_temperature
       endif
    endif

    ! report input parameters
    if (rank == 0) then
       write(30,'(A)') & 
            'Problem: cart-densityfield_fluct (cartesian: density field with fluctuations)'
       write (30,'(A,es10.3)') '1) Average density  (cm^-3): ',average_density_cm3
       write (30,'(2A)') '2) Density fluctuations file (cm^-3): ', &
            densityfluctuations_file
       write (30,'(A)') '3) Temperature setup:'
       if (temperature_setup == 1) then
          write (30,'(A)') 'isobaric'
          write (30,'(A,es10.3)') '4) Average temperature: ',average_temperature
       else
          write (30,'(A)') '  isothermal'
          write (30,'(A,es10.3)') '4) Temperature: ',average_temperature
       endif
    endif
    
#ifdef MPI       
    ! Distribute the input parameters to the other nodes
    ! (densityfield_file is distributed in the MPI ring below)
    call MPI_BCAST(average_density_cm3,1,MPI_DOUBLE_PRECISION,0, & 
         MPI_COMM_NEW,ierror)
    call MPI_BCAST(temperature_setup,1,MPI_INTEGER,0, & 
         MPI_COMM_NEW,ierror)
    call MPI_BCAST(average_temperature,1,MPI_DOUBLE_PRECISION,0, & 
         MPI_COMM_NEW,ierror)
#endif
    
    ! Read in density field (MPI ring)
    if (rank.eq.0) then
       
       open(unit=40,file=densityfluctuations_file,form='UNFORMATTED',status='OLD')
       read(40) mx1,my1,mz1
       write(30,*)  mx1,my1,mz1
       if (mx1 .ne. meshx .or. my1.ne.meshy .or. mz1.ne.meshz) then
          write(30,*) 'Error: ', &
               'density field does not match specified mesh size'
          stop
       endif
       allocate(tempdens(1-mbc:mx1+mbc,1-mbc:my1+mbc,1-mbc:mz1+mbc))
       read(40) tempdens(1:mx1,1:my1,1:mz1)
       close(40)
       
#ifdef MPI
       ! Send filename to next processor
       if (npr > 1) then
          call MPI_ISSEND(densityfluctuations_file,512,MPI_CHARACTER,rank+1, &
               outputcircle,MPI_COMM_NEW,request,ierror)
          write(30,*) rank,'sent filename to ',rank+1
          ! Wait for the circle to complete
          call MPI_RECV(densityfluctuations_file,512,MPI_CHARACTER,npr-1, &
               outputcircle,MPI_COMM_NEW,status,ierror)
          write(30,*) rank,'received filename from ',npr-1
          flush(30)
       endif

    else
       
       ! Receive filename from previous processor
       call MPI_RECV(densityfluctuations_file,512,MPI_CHARACTER,rank-1, &
            outputcircle,MPI_COMM_NEW,status,ierror)
       write(30,*) rank,'received filename from ',rank-1
       flush(30)
       
       if (ierror == 0) then ! if ok
          open(unit=40,file=densityfluctuations_file,form='UNFORMATTED', &
               status='OLD')
          read(40) mx1,my1,mz1
          write(30,*)  mx1,my1,mz1
          if (mx1 .ne. meshx .or. my1.ne.meshy .or. mz1.ne.meshz) then
             write(30,*) 'Error: ', &
                  'density field does not match specified mesh size'
             stop
          endif
          allocate(tempdens(1-mbc:mx1+mbc,1-mbc:my1+mbc,1-mbc:mz1+mbc))
          read(40) tempdens(1:mx1,1:my1,1:mz1)
          close(40)
          
          ! Determine next processor (npr -> 0)
          nextproc=mod(rank+1,npr)
          ! Send filename along
          call MPI_ISSEND(densityfluctuations_file,512,MPI_CHARACTER,nextproc, &
               outputcircle,MPI_COMM_NEW,request,ierror)
          write(30,*) rank,'sent filename to ',nextproc
       else
          write(30,*) 'MPI error in output routine'
          ierror=ierror+1
       endif
#endif
    endif

    flush(30)
    ! Assume outflow boundaries
    do k=1-mbc,mz1+mbc
       do j=1-mbc,my1+mbc
          do i=1-mbc,0
             tempdens(i,j,k)=tempdens(1,j,k)
          enddo
       enddo
    enddo
    do k=1-mbc,mz1+mbc
       do j=1-mbc,0
          do i=1-mbc,mx1+mbc
             tempdens(i,j,k)=tempdens(i,1,k)
          enddo
       enddo
    enddo
    do k=1-mbc,0
       do j=1-mbc,my1+mbc
          do i=1-mbc,mx1+mbc
             tempdens(i,j,k)=tempdens(i,j,1)
          enddo
       enddo
    enddo
    
    ! Set lower fluctuation floor
    do k=1-mbc,mz1+mbc
       do j=1-mbc,my1+mbc
          do i=1-mbc,mx1+mbc
             if (1.0+tempdens(i,j,k) < minimum_fluctuation) &
                  tempdens(i,j,k)=minimum_fluctuation-1.0
          enddo
       enddo
    enddo

    ! Find maximum density (for isobaric initial conditions)
    write(30,*) maxval(tempdens),minval(tempdens),scdens
    flush(30)
    maxdens=maxval(tempdens)!*mu*m_p/scdens 
    
    ! Set density in state array
    ! Scale to program scaling
    do k=sz-mbc,ez+mbc
       do j=sy-mbc,ey+mbc
          do i=sx-mbc,ex+mbc
             ! densities are in cm^-3
             state(i,j,k,RHO)=n2rho(average_density_cm3*(1.0+real(tempdens(i,j,k),dp)))/scdens 
          enddo
       enddo
    enddo
    deallocate(tempdens)
    
    ! Set the initial momenta
    do k=sz-mbc,ez+mbc
       do j=sy-mbc,ey+mbc
          do i=sx-mbc,ex+mbc
             state(i,j,k,RHVX)=0.0d0
             state(i,j,k,RHVY)=0.0d0
             state(i,j,k,RHVZ)=0.0d0
          enddo
       enddo
    enddo
    
    write(30,*) 'setting pressure'
    write(30,*) minval(state(:,:,:,RHO))
    flush(30)
    ! Set the initial pressure and energy density
    ! (depends on the temperature setup, isobaric or isothermal)
    if (temperature_setup == 1) then
       ! isobaric
       epressr=temper2pressr(average_temperature,average_density_cm3,0.0_dp)/SCENER
       do k=sz-mbc,ez+mbc
          do j=sy-mbc,ey+mbc
             do i=sx-mbc,ex+mbc
                pressr(i,j,k)=epressr
                state(i,j,k,EN)=pressr(i,j,k)/gamma1+ &
                     0.5d0*(state(i,j,k,RHVX)*state(i,j,k,RHVX)+ &
                     state(i,j,k,RHVY)*state(i,j,k,RHVY)+ &
                     state(i,j,k,RHVZ)*state(i,j,k,RHVZ))/state(i,j,k,RHO)
                !state(i,j,k,TRACER1)=1.0d0
             enddo
          enddo
       enddo
    else
       ! isothermal
       do k=sz-mbc,ez+mbc
          do j=sy-mbc,ey+mbc
             do i=sx-mbc,ex+mbc
                pressr(i,j,k)=temper2pressr(average_temperature, &
                     rho2n(state(i,j,k,RHO)*scdens),0.0_dp)/SCENER
                state(i,j,k,EN)=pressr(i,j,k)/gamma1+ &
                     0.5d0*(state(i,j,k,RHVX)*state(i,j,k,RHVX)+ &
                     state(i,j,k,RHVY)*state(i,j,k,RHVY)+ &
                     state(i,j,k,RHVZ)*state(i,j,k,RHVZ))/state(i,j,k,RHO)
                !state(i,j,k,TRACER1)=1.0d0
             enddo
          enddo
       enddo
    endif
    
  end subroutine fresh_start_state
  
  !==========================================================================

  subroutine problemboundary (boundary_id,newold)
    
    ! This routine resets the inner boundary to the inflow condition

    ! Version: dummy routine

    integer,intent(in) :: boundary_id
    integer,intent(in) :: newold
    
  end subroutine problemboundary
  
  !==========================================================================

  subroutine apply_grav_force(dt,newold)

    ! Dummy routine

    real(kind=dp),intent(in) :: dt
    integer,intent(in) :: newold

  end subroutine apply_grav_force

end module problem