[FLASH-USERS] Can't find ignored parameter
Aaron Froese
aaron.froese at generalfusion.com
Tue Dec 21 16:17:00 EST 2010
Hi All,
I am trying to modify the Field Loop example problem to have two fluid species, one in each half of the domain. (Currently they have identical properties because I intend to assign a different resistivity to each.) However, I am getting a puzzling error in which the EOS gamma value is first ignored and then the simulation aborts because it can't find it.
RuntimeParameters_read: ignoring unknown parameter "gamma"...
DRIVER_ABORT: ERROR: cannot locate real runtime parameter.
WARNING: requested real parameter 'gamma' not found
I am including the log and Simulation files. Any assistance would be much appreciated.
Best regards,
Aaron Froese
############################################################
############################################################
# log file
############################################################
############################################################
FLASH log file: 12-21-2010 12:32:06.283 Run number: 1
==============================================================================
Number of processors: 8
Dimensionality: 2
Max Number of Blocks/Proc: 1000
Number x zones: 8
Number y zones: 8
Number z zones: 1
Setup stamp: Tue Dec 21 12:29:56 2010
Build stamp: Tue Dec 21 12:31:48 2010
System info: Linux deepthought 2.6.18-164.11.1.el5.centos.plus #1 SMP Wed Jan 20 18:49:35 EST
Version: FLASH 3.3_release
Build directory: /home/aaron/programs/flash3.3/TwoChamberFieldLoop
Setup syntax: /home/aaron/programs/flash3.3/bin/setup.py magnetoHD/TwoChamberFieldLoop -auto -2d +usm -objdir TwoChamberFieldLoop
f compiler flags: /shared/openmpi44/bin/mpif90 -c -O2 -fdefault-real-8 -fdefault-double-8 -ffree-line-length-none -DMAXBLOCKS=1000 -DNXB=8 -DNYB=8 -DNZB=1 -DN_DIM=2
c compiler flags: /shared/openmpi44/bin/mpicc -I /home/aaron/libraries/hdf5-openmpi44/include -DH5_USE_16_API -O2 -c -DMAXBLOCKS=1000 -DNXB=8 -DNYB=8 -DNZB=1 -DN_DIM=2
==============================================================================
Comment: Two Chamber Field Loop MHD problem
==============================================================================
FLASH Units used:
Driver/DriverMain/Unsplit
Grid/GridBoundaryConditions
Grid/GridMain/paramesh/interpolation/Paramesh4/prolong
Grid/GridMain/paramesh/interpolation/prolong
Grid/GridMain/paramesh/paramesh4/Paramesh4.0/PM4_package/headers
Grid/GridMain/paramesh/paramesh4/Paramesh4.0/PM4_package/mpi_source
Grid/GridMain/paramesh/paramesh4/Paramesh4.0/PM4_package/source
Grid/GridMain/paramesh/paramesh4/Paramesh4.0/PM4_package/utilities/multigrid
Grid/localAPI
IO/IOMain/hdf5/serial/PM
IO/localAPI
Multispecies/MultispeciesMain
Particles
PhysicalConstants/PhysicalConstantsMain
RuntimeParameters/RuntimeParametersMain
Simulation/SimulationMain/magnetoHD/TwoChamberFieldLoop
flashUtilities/contiguousConversion
flashUtilities/general
flashUtilities/interpolation/oneDim
flashUtilities/nameValueLL
monitors/Logfile/LogfileMain
monitors/Profiler
monitors/Timers/TimersMain/MPINative
physics/Cosmology
physics/Diffuse
physics/Eos/EosMain/Multigamma
physics/Gravity
physics/Hydro/HydroMain/unsplit/MHD_StaggeredMesh
physics/materialProperties/Conductivity/ConductivityMain/Constant-diff
physics/materialProperties/MagneticResistivity/MagneticResistivityMain/Constant
physics/materialProperties/MassDiffusivity
physics/materialProperties/Viscosity/ViscosityMain/Constant
physics/sourceTerms/Burn
physics/sourceTerms/Cool
physics/sourceTerms/Heat
physics/sourceTerms/Ionize
physics/sourceTerms/Stir
==============================================================================
RuntimeParameters:
==============================================================================
bndpriorityone = 1
bndprioritythree = 3
bndprioritytwo = 2
checkpointfileintervalstep = 0
checkpointfilenumber = 0
dr_abortpause = 2
fileformatversion = 9
forcedplotfilenumber = 0
gr_restrictallmethod = 3
interpol_order = 2
irenorm = 0
lrefine_del = 0
lrefine_max = 6 [CHANGED]
lrefine_max_prev = 1
lrefine_min = 1
max_particles_per_blk = 100
memory_stat_freq = 100000
min_particles_per_blk = 1
nbegin = 1
nblockx = 1
nblocky = 1
nblockz = 1
nend = 1000000 [CHANGED]
nrefs = 2
order = 2
outputsplitnum = 1
plotfileintervalstep = 0
plotfilenumber = 0
refine_var_count = 4
rolling_checkpoint = 10000
sweeporder = 123
transorder = 3
visc_whichcoefficientisconst = 2
wr_integrals_freq = 1
az_initial = 0.100E-02
limitedslopebeta = 0.100E+01
r_fieldloop = 0.300E+00
cfl = 0.800E+00
checkpointfileintervaltime = 0.100E+00 [CHANGED]
checkpointfileintervalz = 0.180+309
derefine_cutoff_1 = 0.200E+00
derefine_cutoff_2 = 0.200E+00
derefine_cutoff_3 = 0.200E+00
derefine_cutoff_4 = 0.200E+00
diff_constant = 0.100E-01 [CHANGED]
diff_visc_mu = 0.100E+00
diff_visc_nu = 0.400E-03 [CHANGED]
dtinit = 0.100E-09
dtmax = 0.100E+06
dtmin = 0.100E-09
eintswitch = 0.100E-05 [CHANGED]
gr_lrefinemaxredlogbase = 0.100E+02
gr_lrefinemaxredradiusfact = 0.000E+00
gr_lrefinemaxredtref = 0.000E+00
gr_lrefinemaxredtimescale = 0.100E+01
plotfileintervaltime = 0.100E+01
plotfileintervalz = 0.180+309
refine_cutoff_1 = 0.800E+00
refine_cutoff_2 = 0.800E+00
refine_cutoff_3 = 0.800E+00
refine_cutoff_4 = 0.800E+00
refine_filter_1 = 0.100E-01
refine_filter_2 = 0.100E-01
refine_filter_3 = 0.100E-01
refine_filter_4 = 0.100E-01
resistivity = 0.100E-03 [CHANGED]
sim_divide = 0.500E+00
sim_pleft = 0.100E+01
sim_pright = 0.100E+01
sim_rholeft = 0.100E+01
sim_rhoright = 0.000E+00
sim_uleft = 0.100E+01
sim_uright = 0.100E+01
small = 0.100E-09
smalle = 0.100E-09
smallp = 0.100E-09
smallt = 0.100E-09
smallu = 0.100E-09
smallx = 0.100E-09
smlrho = 0.100E-09
tinitial = 0.000E+00
tiny = 0.100E-15
tmax = 0.200E+01 [CHANGED]
tstep_change_factor = 0.200E+01
viscsuppressfactor = 0.100E+01
visctemphigh = 0.150E+09
visctemplow = 0.300E+06
wall_clock_checkpoint = 0.432E+05
wall_clock_time_limit = 0.605E+06
xctr = 0.500E+00
x_refine_center = 0.000E+00
xmax = 0.100E+01
xmin = 0.000E+00
yctr = 0.500E+00
y_refine_center = 0.000E+00
ymax = 0.200E+01 [CHANGED]
ymin = 0.000E+00
zctr = 0.000E+00
zfinal = 0.000E+00
zinitial = -0.100E+01
z_refine_center = 0.000E+00
zmax = 0.100E+01
zmin = 0.000E+00
riemannsolver = Roe
unitsystem = none
basenm = twoChamberFieldLoop_ [CHANGED]
entropyfixmethod = HARTENHYMAN
eosmode = dens_ie
eosmodeinit = dens_ie
geometry = cartesian
grav_boundary_type = isolated
log_file = twoChamberFieldLoop.log [CHANGED]
output_directory =
pc_unitsbase = CGS
plot_grid_var_1 = none
plot_grid_var_10 = none
plot_grid_var_11 = none
plot_grid_var_12 = none
plot_grid_var_2 = none
plot_grid_var_3 = none
plot_grid_var_4 = none
plot_grid_var_5 = none
plot_grid_var_6 = none
plot_grid_var_7 = none
plot_grid_var_8 = none
plot_grid_var_9 = none
plot_var_1 = DEUT dens [CHANGED]
plot_var_10 = magp [CHANGED]
plot_var_11 = none
plot_var_12 = none
plot_var_13 = none
plot_var_14 = none
plot_var_15 = none
plot_var_16 = none
plot_var_17 = none
plot_var_18 = none
plot_var_19 = none
plot_var_2 = DEUT pres [CHANGED]
plot_var_3 = DEUT temp [CHANGED]
plot_var_4 = TRIT dens [CHANGED]
plot_var_5 = TRIT pres [CHANGED]
plot_var_6 = TRIT temp [CHANGED]
plot_var_7 = divb [CHANGED]
plot_var_8 = magx [CHANGED]
plot_var_9 = magy [CHANGED]
prof_file = profile.dat
prolmethod = injection_prol [CHANGED]
refine_var_1 = magx [CHANGED]
refine_var_2 = magy [CHANGED]
refine_var_3 = none
refine_var_4 = none
run_comment = Two Chamber Field Loop MHD pro [CHANGED]
run_number = 1
slopelimiter = mc [CHANGED]
stats_file = flash.dat
xl_boundary_type = periodic
xr_boundary_type = periodic
yl_boundary_type = periodic
yr_boundary_type = periodic
zl_boundary_type = periodic
zr_boundary_type = periodic
eosforriemann = F
e_modification = T
forcehydrolimit = F
alwayscomputeuservars = T
alwaysrestrictcheckpoint = T
bytepack = F
charlimiting = T
chkguardcellsinput = F
chkguardcellsoutput = F
converttoconsvdformeshcalls = F
converttoconsvdinmeshinterp = T
corners = F
dr_shortenlaststepbeforetma = F
eachprocwritesownabortlog = F
eachprocwritessummary = T
earlyblockdistadjustment = T
enablemaskedgcfill = F
energyfix = T [CHANGED]
entropy = F
facevar2ndorder = T
flux_correct = T
fullylimit = F
geometryoverride = F
gr_lrefinemaxreddobylogr = F
gr_lrefinemaxreddobytime = F
ignoreforcedplot = F
killdivb = T
plotfilegridquantitydp = F
plotfilemetadatadp = F
refine_on_particle_count = F
restart = F
shockdetect = F
shockinstabilityfix = F
unbiased_geometry = F
useburn = F
usecollectivehdf5 = F
useconductivity = F [CHANGED]
usecool = F
usediffuse = F
usegravity = F
useheat = F
usehydro = T
uselegacylabels = T
usemagneticresistivity = F [CHANGED]
usemassdiffusivity = F
useparticles = F
usestir = T
useviscosity = F [CHANGED]
use_cma_advection = F
use_cma_flattening = F
use_cma_steepening = F
use_flattening = F
use_gravconsv = F
use_gravhalfupdate = F
use_steepening = F
writestatsummary = T
==============================================================================
WARNING: Ignored Parameters :
These parameters were found in the flash.par file, but they were
not declared in any Config file for the simulation!
gamma
==============================================================================
Known units of measurement:
Unit CGS Value Base Unit
1 cm 1.0000 cm
2 s 1.0000 s
3 K 1.0000 K
4 g 1.0000 g
5 esu 1.0000 esu
6 m 100.00 cm
7 km 0.10000E+06 cm
8 pc 0.30857E+19 cm
9 kpc 0.30857E+22 cm
10 Mpc 0.30857E+25 cm
11 Gpc 0.30857E+28 cm
12 Rsun 0.69600E+11 cm
13 AU 0.14960E+14 cm
14 yr 0.31557E+08 s
15 Myr 0.31557E+14 s
16 Gyr 0.31557E+17 s
17 kg 1000.0 g
18 Msun 0.19889E+34 g
19 amu 0.16605E-23 g
20 C 0.29979E+10 esu
21 LFLY 0.30857E+25 cm
22 TFLY 0.20576E+18 s
23 MFLY 0.98847E+46 g
24 clLength 0.30857E+25 cm
25 clTime 0.31557E+17 s
26 clMass 0.19889E+49 g
27 clTemp 0.11604E+08 K
-----------End of Units--------------------
Known physical constants:
Constant Name Constant Value cm s g K esu
1 Newton 0.66726E-07 3.00 -2.00 -1.00 0.00 0.00
2 speed of light 0.29979E+11 1.00 -1.00 0.00 0.00 0.00
3 Planck 0.66261E-26 2.00 -1.00 1.00 0.00 0.00
4 electron charge 0.48032E-09 0.00 0.00 0.00 0.00 1.00
5 electron mass 0.91094E-27 0.00 0.00 1.00 0.00 0.00
6 proton mass 0.16726E-23 0.00 0.00 1.00 0.00 0.00
7 fine-structure 0.72974E-02 0.00 0.00 0.00 0.00 0.00
8 Avogadro 0.60221E+24 0.00 0.00 0.00 0.00 0.00
9 Boltzmann 0.13807E-15 2.00 -2.00 1.00 -1.00 0.00
10 ideal gas constant 0.83145E+08 2.00 -2.00 0.00 -1.00 0.00
11 Wien 0.28978 1.00 0.00 0.00 1.00 0.00
12 Stefan-Boltzmann 0.56705E-04 0.00 -3.00 1.00 -4.00 0.00
13 pi 3.1416 0.00 0.00 0.00 0.00 0.00
14 e 2.7183 0.00 0.00 0.00 0.00 0.00
15 Euler 0.57722 0.00 0.00 0.00 0.00 0.00
==============================================================================
Multifluid database contents:
Initially defined values of species:
Name Index Total Positive Neutral Negative bind Ener Gamma
deut 18 2.00E+00 1.00E+00 -9.99E+02 -9.99E+02 -9.99E+02 1.67E+00
trit 19 3.00E+00 1.00E+00 -9.99E+02 -9.99E+02 -9.99E+02 1.67E+00
==============================================================================
[ 12-21-2010 12:32:06.384 ] [gr_initGeometry] checking BCs for idir: 1
WARNING: requested real parameter 'gamma' not found
[ 12-21-2010 12:32:06.400 ] [gr_initGeometry] checking BCs for idir: 2
WARNING: requested real parameter 'gamma' not found
WARNING: requested real parameter 'gamma' not found
WARNING: requested real parameter 'gamma' not found
WARNING: requested real parameter 'gamma' not found
[ 12-21-2010 12:32:06.457 ] [DRIVER_ABORT]: Driver_abort() called by PE 0
WARNING: requested real parameter 'gamma' not found
[ 12-21-2010 12:32:06.483 ] abort_message: ERROR: cannot locate real runtime parameter.
WARNING: requested real parameter 'gamma' not found
WARNING: requested real parameter 'gamma' not found
############################################################
############################################################
# Config
############################################################
############################################################
# Configuration file for 2D Field Loop MHD problem
# James Stone, www.astro.princeton.edu/~jstone
REQUIRES Driver
REQUIRES Grid
REQUESTS IO
REQUIRES Multispecies/MultispeciesMain
REQUIRES physics/Hydro/HydroMain
REQUIRES physics/Eos/EosMain/Multigamma
# ----------------- For Resistive MHD setup ------------------#
REQUIRES physics/materialProperties/Conductivity/ConductivityMain/Constant-diff
REQUIRES physics/materialProperties/Viscosity/ViscosityMain
REQUIRES physics/materialProperties/MagneticResistivity/MagneticResistivityMain
#------------------ End of Resistive MHD setup ----------------#
GUARDCELLS 4
D tiny Threshold value used for numerical zero
PARAMETER tiny REAL 1.e-16
D xCtr x center of field loop
PARAMETER xCtr REAL 0.5
D yCtr y center of field loop
PARAMETER yCtr REAL 0.5
D zCtr z center of field loop
PARAMETER zCtr REAL 0.0
D Az_initial Initial z-component of magnetic vector potential
PARAMETER Az_initial REAL 0.001
D R_fieldLoop Radius of field loop
PARAMETER R_fieldLoop REAL 0.3
D sim_rhoLeft Density in the left part of the grid
PARAMETER sim_rhoLeft REAL 1. [0 to 1]
D sim_rhoRight Density in the right part of the grid
PARAMETER sim_rhoRight REAL 0. [0 to 1]
D sim_pLeft Pressure in the left part of the grid
PARAMETER sim_pLeft REAL 1. [0 to ]
D sim_pRight Pressure in the righ part of the grid
PARAMETER sim_pRight REAL 1. [0 to ]
D sim_uLeft fluid velocity in the left part of the grid
PARAMETER sim_uLeft REAL 1.
D sim_uRight fluid velocity in the right part of the grid
PARAMETER sim_uRight REAL 1.
D sim_divide Point of intersection between the shock plane and the x-axis
PARAMETER sim_divide REAL 0.5
VARIABLE vecz # vector potential Az
VARIABLE curz # vector potential jz
SPECIES DEUT
SPECIES TRIT
############################################################
############################################################
# flash.par
############################################################
############################################################
# Runtime parameters for the MHD ResistiveFieldLoop problem.
# Please also look at default.par for more runtime parameters in the object directory
# Grid dimensionality and geometry
geometry = "cartesian"
# Size of computational volume for 2D test
xmin = 0.
xmax = 1.
ymin = 0.
ymax = 2.
xCtr = 0.5
yCtr = 0.5
# Size of computational volume for 3D test
#xmin =-0.5
#xmax = 0.5
#ymin =-0.5
#ymax = 0.5
#zmin =-1.0
#zmax = 1.0
# Density, pressure, and velocity on either side of interface
sim_rhoLeft = 1.0
sim_rhoRight = 0.0
sim_pLeft = 1.0
sim_pRight = 1.0
sim_uLeft = 1.0
sim_uRight = 1.0
sim_divide = 0.5
# Gas ratio of specific heats
gamma = 1.666666667
# Initial strength of the magnetic vector potential
Az_initial = 0.001
# Radius of the field loop
R_fieldLoop = 0.3
# Boundary conditions (code -22 is periodic)
xl_boundary_type = "periodic"
xr_boundary_type = "periodic"
yl_boundary_type = "periodic"
yr_boundary_type = "periodic"
zl_boundary_type = "periodic"
zr_boundary_type = "periodic"
# Simulation (grid, time, I/O) parameters
run_comment = "Two Chamber Field Loop MHD problem"
log_file = "twoChamberFieldLoop.log"
basenm = "twoChamberFieldLoop_"
restart = .false.
#checkPointFileNumber=1
#plotFileNumber = 1
nend = 1000000
tmax = 2.0
cfl = 0.8
plot_var_1 = "DEUT dens"
plot_var_2 = "DEUT pres"
plot_var_3 = "DEUT temp"
plot_var_4 = "TRIT dens"
plot_var_5 = "TRIT pres"
plot_var_6 = "TRIT temp"
plot_var_7 = "divb"
plot_var_8 = "magx"
plot_var_9 = "magy"
plot_var_10 = "magp"
convertToConsvdInMeshInterp = .true.
checkpointFileIntervalTime = 0.1
#checkpointFileIntervalStep = 10
# AMR parameters
#nblockx = 1
#nblocky = 1
lrefine_min = 1
lrefine_max = 6
nrefs = 2
refine_var_1 = "magx"
refine_var_2 = "magy"
eintSwitch = 1.e-6
# DivB control switch
killdivb = .true.
# Flux Conservation for AMR
flux_correct = .true.
## -------------------------------------------------------------##
## SWITCHES SPECIFIC TO THE UNSPLIT STAGGERED MESH MHD SOLVER ##
# I. INTERPOLATION SCHEME:
order = 2 # Interpolation order (First/Second order)
slopeLimiter = "mc" # Slope limiters (minmod, mc, vanLeer, hybrid, limited)
LimitedSlopeBeta= 1. # Slope parameter for the "limited" slope by Toro
charLimiting = .true. # Characteristic limiting vs. Primitive limiting
# II. MAGNETIC(B) AND ELECTRIC(E) FIELDS:
E_modification = .true. # High order algorithm for E-field construction
energyFix = .true. # Update magnetic energy using staggered B-fields
ForceHydroLimit = .false. # Pure Hydro Limit (B=0)
prolMethod = "injection_prol" # Prolongation method (injecton_prol, balsara_prol)
# III. RIEMANN SOLVERS:
RiemannSolver = "Roe" # LLF, HLL, HLLC, HLLD, Marquina, Roe
shockInstabilityFix = .false. # Carbuncle instability fix for the Roe solver
entropy = .false. # Entropy fix for the Roe solver
# IV. STRONG SHOCK HANDELING SCHEME:
shockDetect = .false. # Shock Detect for numerical stability
## -------------------------------------------------------------##
## -------------------------------------------------------##
## Resistive MHD parameters ##
useMagneticResistivity = .false.
resistivity = 1.0E-4 # magnetic Resistivity
useViscosity = .false.
diff_visc_nu = 4.0E-4
useConductivity = .false.
diff_constant = 1.0E-2
## -------------------------------------------------------##
############################################################
############################################################
# Simulation_data.F90
############################################################
############################################################
module Simulation_data
implicit none
#include "constants.h"
!! *** Runtime Parameters *** !!
real, save :: sim_gamma
real, save :: sim_smallP
real, save :: sim_smallX
real, save :: sim_xMin, sim_xMax, sim_yMin, sim_yMax, sim_zMin, sim_zMax
real, save :: sim_xCtr, sim_yCtr, sim_zCtr, sim_divide
real, save :: sim_Az_initial, sim_fieldLoopRadius
real, save :: sim_rhoLeft, sim_rhoRight
real, save :: sim_pLeft, sim_pRight
real, save :: sim_uLeft, sim_uRight
logical, save :: sim_gCell, sim_killdivb
end module Simulation_data
############################################################
############################################################
# Simulation_initBlock.F90
############################################################
############################################################
subroutine Simulation_initBlock(blockID, myPE)
!use Simulation_speciesData, ONLY : sim_specElementSymbol,sim_specElement,&
! sim_specSelected, sim_specNumElect, SPEC_NUM, SPEC_NIMAX, sim_specXfrac
use Simulation_data, ONLY : sim_gCell,sim_xCtr,sim_yCtr,sim_fieldLoopRadius, &
& sim_Az_initial,sim_smallP,sim_gamma &
& sim_divide,sim_rhoLeft,sim_pLeft,sim_uLeft, &
& sim_rhoRight,sim_pRight,sim_uRight,sim_smallX
use Grid_interface, ONLY : Grid_getBlkIndexLimits, &
Grid_getCellCoords, &
Grid_getDeltas, &
Grid_getBlkPtr, &
Grid_releaseBlkPtr
use Simulation_data, ONLY : sim_killdivb
implicit none
#include "constants.h"
#include "Flash.h"
#include "Eos.h"
#include "Multispecies.h"
!! Arguments ------------------------
integer, intent(in) :: blockID, myPE
!! ----------------------------------
integer :: i, j, k, n, istat, sizeX, sizeY, sizeZ
integer, dimension(2,MDIM) :: blkLimits, blkLimitsGC
real :: rhoZone, velxZone, velyZone, velzZone, presZone, &
eintZone, enerZone, ekinZone
real :: radius, dx, dy, dz
real, allocatable,dimension(:) :: xCoord,xCoordL,xCoordR,&
yCoord,yCoordL,yCoordR,&
zCoord,zCoordL,zCoordR
real, dimension(MDIM) :: del
real, pointer, dimension(:,:,:,:) :: solnData, facexData,faceyData,facezData
real :: xx,xxL,xxR,yy,zz
#ifdef FIXEDBLOCKSIZE
real, dimension(GRID_IHI_GC+1,GRID_JHI_GC+1,GRID_KHI_GC+1) :: Az,Ax,Ay
#else
real, allocatable, dimension(:,:,:) :: Az,Ax,Ay
#endif
! dump some output to stdout listing the paramters
if (myPE == MASTER_PE) then
1 format (1X, 1P, 4(A7, E13.7, :, 1X))
2 format (1X, 1P, 2(A7, E13.7, 1X), A7, I13)
endif
! get the integer index information for the current block
call Grid_getBlkIndexLimits(blockId,blkLimits,blkLimitsGC)
sizeX = blkLimitsGC(HIGH,IAXIS)-blkLimitsGC(LOW,IAXIS)+1
sizeY = blkLimitsGC(HIGH,JAXIS)-blkLimitsGC(LOW,JAXIS)+1
sizeZ = blkLimitsGC(HIGH,KAXIS)-blkLimitsGC(LOW,KAXIS)+1
allocate(xCoord(sizeX), stat=istat)
allocate(xCoordL(sizeX),stat=istat)
allocate(xCoordR(sizeX),stat=istat)
allocate(yCoord(sizeY), stat=istat)
allocate(yCoordL(sizeY),stat=istat)
allocate(yCoordR(sizeY),stat=istat)
allocate(zCoord(sizeZ), stat=istat)
allocate(zCoordL(sizeZ),stat=istat)
allocate(zCoordR(sizeZ),stat=istat)
xCoord = 0.0
xCoordL = 0.0
xCoordR = 0.0
yCoord = 0.0
yCoordL = 0.0
yCoordR = 0.0
zCoord = 0.0
zCoordL = 0.0
zCoordR = 0.0
#ifndef FIXEDBLOCKSIZE
if (NDIM == 2) then
allocate(Ax(sizeX+1,sizeY+1,1),stat=istat)
allocate(Ay(sizeX+1,sizeY+1,1),stat=istat)
allocate(Az(sizeX+1,sizeY+1,1),stat=istat)
elseif (NDIM == 3) then
allocate(Ax(sizeX+1,sizeY+1,sizeZ+1),stat=istat)
allocate(Ay(sizeX+1,sizeY+1,sizeZ+1),stat=istat)
allocate(Az(sizeX+1,sizeY+1,sizeZ+1),stat=istat)
endif
#endif
if (NDIM == 3) then
call Grid_getCellCoords(KAXIS,blockId,CENTER, sim_gCell,zCoord, sizeZ)
call Grid_getCellCoords(KAXIS,blockId,LEFT_EDGE, sim_gCell,zCoordL,sizeZ)
call Grid_getCellCoords(KAXIS,blockId,RIGHT_EDGE,sim_gCell,zCoordR,sizeZ)
endif
if (NDIM >= 2) then
call Grid_getCellCoords(JAXIS,blockId,CENTER, sim_gCell,yCoord, sizeY)
call Grid_getCellCoords(JAXIS,blockId,LEFT_EDGE, sim_gCell,yCoordL,sizeY)
call Grid_getCellCoords(JAXIS,blockId,RIGHT_EDGE,sim_gCell,yCoordR,sizeY)
endif
call Grid_getCellCoords(IAXIS,blockId,CENTER, sim_gCell,xCoord, sizeX)
call Grid_getCellCoords(IAXIS,blockId,LEFT_EDGE, sim_gCell,xCoordL,sizeX)
call Grid_getCellCoords(IAXIS,blockId,RIGHT_EDGE,sim_gCell,xCoordR,sizeX)
!------------------------------------------------------------------------------
! Construct Az at each cell corner
! Bx = dAz/dy - dAy/dz
! By = dAx/dz - dAz/dx
! Bz = dAy/dx - dAx/dy
Az = 0.
Ax = 0.
Ay = 0.
call Grid_getDeltas(blockID,del)
dx = del(1)
dy = del(2)
dz = del(3)
if (NDIM == 2) then
k = 1
do j = blkLimitsGC(LOW,JAXIS),blkLimitsGC(HIGH,JAXIS)+1
do i = blkLimitsGC(LOW,IAXIS),blkLimitsGC(HIGH,IAXIS)+1
#if NFACE_VARS > 1
! x Coord at cell corner
if (i <=blkLimitsGC(HIGH,IAXIS)) then
xx = xCoordL(i)
else
xx = xCoordR(i-1)
endif
! y Coord at cell corner
if (j <=blkLimitsGC(HIGH,JAXIS)) then
yy = yCoordL(j)
else
yy = yCoordR(j-1)
endif
#else
! x Coord at cell center
if (i <=blkLimitsGC(HIGH,IAXIS)) then
xx = xCoord(i)
else
xx = xCoord(i-1) + dx
endif
! y Coord at cell center
if (j <=blkLimitsGC(HIGH,JAXIS)) then
yy = yCoord(j)
else
yy = yCoord(j-1) + dy
endif
#endif
! define radius of the field loop
radius = sqrt((xx-sim_xCtr)**2 + (yy-sim_yCtr)**2)
if (radius <= sim_fieldLoopRadius) then
Ax(i,j,k) = 0.
Ay(i,j,k) = 0.
Az(i,j,k) = sim_Az_initial*(sim_fieldLoopRadius - radius)
else
Ax(i,j,k) = 0.
Ay(i,j,k) = 0.
Az(i,j,k) = 0.
endif
enddo
enddo
elseif (NDIM == 3) then
do k = blkLimitsGC(LOW,KAXIS),blkLimitsGC(HIGH,KAXIS)+1
do j = blkLimitsGC(LOW,JAXIS),blkLimitsGC(HIGH,JAXIS)+1
do i = blkLimitsGC(LOW,IAXIS),blkLimitsGC(HIGH,IAXIS)+1
!! Rotated coordinate system
!! / xx \ / cos 0 sin \
!! | yy | = | 0 1 0 |
!! \ zz / \ -sin 0 cos /
#if NFACE_VARS > 0
! x Coord at cell corner
if (i <=blkLimitsGC(HIGH,IAXIS)) then
xx = xCoordL(i)
else
xx = xCoordR(i-1)
endif
! y Coord at cell corner
if (j <=blkLimitsGC(HIGH,JAXIS)) then
yy = yCoordL(j)
else
yy = yCoordR(j-1)
endif
! z Coord at cell corner
if (k <=blkLimitsGC(HIGH,KAXIS)) then
zz = zCoordL(k)
else
zz = zCoordR(k-1)
endif
#else
! x Coord at cell center
if (i <=blkLimitsGC(HIGH,IAXIS)) then
xx = xCoord(i)
else
xx = xCoord(i-1) + dx
endif
! y Coord at cell center
if (j <=blkLimitsGC(HIGH,JAXIS)) then
yy = yCoord(j)
else
yy = yCoord(j-1) + dy
endif
! z Coord at cell center
if (k <=blkLimitsGC(HIGH,KAXIS)) then
zz = zCoord(k)
else
zz = zCoord(k-1) + dz
endif
#endif
!xx = (cos_ang*xx + sin_ang*zz) ! with rotation
xx = xx !without any rotation
yy = yy
!do while (xx > 0.5*lambda)
! xx = xx - lambda
!enddo
!do while (xx < -0.5*lambda)
! xx = xx + lambda
!enddo
radius = sqrt(xx**2 + yy**2)
if (radius <= sim_fieldLoopRadius) then
Ax(i,j,k) = 0.
Ay(i,j,k) = 0.
Az(i,j,k) = sim_Az_initial*(sim_fieldLoopRadius - radius)
else
Ax(i,j,k) = 0.
Ay(i,j,k) = 0.
Az(i,j,k) = 0.
endif
enddo
enddo
enddo
endif
call Grid_getBlkPtr(blockID,solnData,CENTER)
#if NFACE_VARS > 0
if (sim_killdivb) then
call Grid_getBlkPtr(blockID,facexData,FACEX)
call Grid_getBlkPtr(blockID,faceyData,FACEY)
if (NDIM == 3) call Grid_getBlkPtr(blockID,facezData,FACEZ)
endif
#endif
!------------------------------------------------------------------------------
! Loop over cells in the block. For each, compute the physical position of
! its left and right edge and its center as well as its physical width.
! Then decide which side of the initial discontinuity it is on and initialize
! the hydro variables appropriately.
do k = blkLimitsGC(LOW,KAXIS),blkLimitsGC(HIGH,KAXIS)
! get the coordinates of the cell center in the z-direction
zz = zCoord(k)
do j = blkLimitsGC(LOW,JAXIS),blkLimitsGC(HIGH,JAXIS)
! get the coordinates of the cell center in the y-direction
yy = yCoord(j)
do i = blkLimitsGC(LOW,IAXIS),blkLimitsGC(HIGH,IAXIS)
! get the cell center, left, and right positions in x
xx = xCoord(i)
xxL = xCoordL(i)
xxR = xCoordR(i)
! initialize cells to the left of the initial shock.
if (xxR <= sim_divide) then
rhoZone = sim_rhoLeft
presZone = sim_pLeft
velxZone = 0.0
velyZone = sim_uLeft
velzZone = 0.0
! initialize cells which straddle the shock. Treat them as though 1/2 of
! the cell lay to the left and 1/2 lay to the right.
elseif ((xxL < sim_divide) .and. (xxR > sim_divide)) then
rhoZone = 0.5 * (sim_rhoLeft+sim_rhoRight)
presZone = 0.5 * (sim_pLeft+sim_pRight)
velxZone = 0.0
velyZone = 0.5 *(sim_uLeft+sim_uRight)
velzZone = 0.0
! initialize cells to the right of the initial shock.
else
rhoZone = sim_rhoRight
presZone = sim_pRight
velxZone = 0.0
velyZone = sim_uRight
velzZone = 0.0
endif
if( rhoZone < sim_smallX ) then
rhoZone = sim_smallX
endif
!assign species fractions
solnData(SPECIES_BEGIN ,i,j,k) = rhoZone
solnData(SPECIES_BEGIN+1,i,j,k) = 1.0e0-rhoZone
! Compute the gas energy and set the gamma-values needed for the equation of
! state.
ekinZone = 0.5 * (velxZone**2 + velyZone**2 + velzZone**2)
! specific internal energy
eintZone = presZone / ((sim_gamma-1.)*rhoZone)
! total specific gas energy
enerZone = eintZone + ekinZone
! Take a limit value
enerZone = max(enerZone, sim_smallP)
! store the variables in the current zone via Grid put methods
! data is put stored one cell at a time with these calls to Grid_putData
solnData(DENS_VAR,i,j,k)= 1.0
solnData(VELX_VAR,i,j,k)= velxZone
solnData(VELY_VAR,i,j,k)= velyZone
solnData(VELZ_VAR,i,j,k)= velzZone
solnData(PRES_VAR,i,j,k)= presZone
solnData(TEMP_VAR,i,j,k)= 1.0
#ifdef ENER_VAR
solnData(ENER_VAR,i,j,k)= enerZone
#endif
#ifdef EINT_VAR
solnData(EINT_VAR,i,j,k)= eintZone
#endif
#ifdef GAME_VAR
solnData(GAME_VAR,i,j,k)= sim_gamma
#endif
#ifdef GAMC_VAR
solnData(GAMC_VAR,i,j,k)= sim_gamma
#endif
enddo
enddo
enddo
!------------------------------------------------------
do k = blkLimitsGC(LOW,KAXIS),blkLimitsGC(HIGH,KAXIS)
do j = blkLimitsGC(LOW,JAXIS),blkLimitsGC(HIGH,JAXIS)
do i = blkLimitsGC(LOW,IAXIS),blkLimitsGC(HIGH,IAXIS)
#ifdef VECZ_VAR
! vector potential Az
if (NFACE_VARS > 1) then
solnData(VECZ_VAR,i,j,k) = .25*(Az(i,j,k)+Az(i+1,j,k)+Az(i,j+1,k)+Az(i+1,j+1,k))
else
solnData(VECZ_VAR,i,j,k) = Az(i,j,k)
endif
#endif
#if NFACE_VARS > 0
!! In this case we initialized Az using the cell-cornered coordinates.
if (sim_killdivb) then
if (NDIM == 2) then
facexData(MAG_FACE_VAR,i,j,k)= (Az(i,j+1,k)-Az(i,j,k))/dy
faceyData(MAG_FACE_VAR,i,j,k)=-(Az(i+1,j,k)-Az(i,j,k))/dx
elseif (NDIM == 3) then
facexData(MAG_FACE_VAR,i,j,k)= -(Ay(i,j,k+1)-Ay(i,j,k))/dz + (Az(i,j+1,k)-Az(i,j,k))/dy
faceyData(MAG_FACE_VAR,i,j,k)= (Ax(i,j,k+1)-Ax(i,j,k))/dz - (Az(i+1,j,k)-Az(i,j,k))/dx
facezData(MAG_FACE_VAR,i,j,k)= -(Ax(i,j+1,k)-Ax(i,j,k))/dy + (Ay(i+1,j,k)-Ay(i,j,k))/dx
endif
endif
#else
!! In this case we initialized Az using the cell-centered coordinates.
if (NDIM == 2) then
solnData(MAGX_VAR,i,j,k)= 0.5*(Az(i,j+1,k)-Az(i,j-1,k))/dy
solnData(MAGY_VAR,i,j,k)=-0.5*(Az(i+1,j,k)-Az(i-1,j,k))/dx
elseif (NDIM == 3) then
solnData(MAGX_VAR,i,j,k)= -0.5*((Ay(i,j,k+1)-Ay(i,j,k-1))/dz + (Az(i,j+1,k)-Az(i,j-1,k))/dy)
solnData(MAGY_VAR,i,j,k)= 0.5*((Ax(i,j,k+1)-Ax(i,j,k-1))/dz - (Az(i+1,j,k)-Az(i-1,j,k))/dx)
solnData(MAGZ_VAR,i,j,k)= -0.5*((Ax(i,j+1,k)-Ax(i,j-1,k))/dy + (Ay(i+1,j,k)-Ay(i-1,j,k))/dx)
endif
#endif
enddo
enddo
enddo
do k=blkLimitsGC(LOW,KAXIS),blkLimitsGC(HIGH,KAXIS)
do j = blkLimitsGC(LOW,JAXIS),blkLimitsGC(HIGH,JAXIS)
do i = blkLimitsGC(LOW,IAXIS),blkLimitsGC(HIGH,IAXIS)
#if NFACE_VARS > 0
solnData(MAGX_VAR,i,j,k) = 0.5*(facexData(MAG_FACE_VAR,i,j,k)+facexData(MAG_FACE_VAR,i+1,j,k))
solnData(MAGY_VAR,i,j,k) = 0.5*(faceyData(MAG_FACE_VAR,i,j,k)+faceyData(MAG_FACE_VAR,i,j+1,k))
if (NDIM == 3) then
solnData(MAGZ_VAR,i,j,k) = 0.5*(facezData(MAG_FACE_VAR,i,j,k)+facezData(MAG_FACE_VAR,i,j,k+1))
endif
#if NDIM == 1
solnData(DIVB_VAR,i,j,k) = 0.
#elif NDIM >= 2
solnData(DIVB_VAR,i,j,k)= &
(facexData(MAG_FACE_VAR,i+1,j, k ) - facexData(MAG_FACE_VAR,i,j,k))/dx &
+ (faceyData(MAG_FACE_VAR,i, j+1,k ) - faceyData(MAG_FACE_VAR,i,j,k))/dy
#if NDIM == 3
solnData(DIVB_VAR,i,j,k)= solnData(DIVB_VAR,i,j,k) &
+ (facezData(MAG_FACE_VAR,i, j, k+1) - facezData(MAG_FACE_VAR,i,j,k))/dz
#endif
#endif
#else !NFACE_VARS == 0
solnData(DIVB_VAR,i,j,k) = 0.
#endif !NFACE_VARS
! Update the magnetic pressure
solnData(MAGP_VAR,i,j,k) = .5*dot_product(solnData(MAGX_VAR:MAGZ_VAR,i,j,k),&
solnData(MAGX_VAR:MAGZ_VAR,i,j,k))
enddo
enddo
enddo
#ifdef CURJ_VAR
k=1
do j = blkLimitsGC(LOW,JAXIS),blkLimitsGC(HIGH,JAXIS)
do i = blkLimitsGC(LOW,IAXIS),blkLimitsGC(HIGH,IAXIS)
solnData(CURJ_VAR,i,j,k) = &
(faceyData(MAG_FACE_VAR,i+1,j,k)-faceyData(MAG_FACE_VAR,i,j,k))/dx &
-(facexData(MAG_FACE_VAR,i,j+1,k)-facexData(MAG_FACE_VAR,i,j,k))/dy
enddo
enddo
#endif
! Release pointer
call Grid_releaseBlkPtr(blockID,solnData,CENTER)
#if NFACE_VARS > 0
if (sim_killdivb) then
call Grid_releaseBlkPtr(blockID,facexData,FACEX)
call Grid_releaseBlkPtr(blockID,faceyData,FACEY)
if (NDIM == 3) call Grid_releaseBlkPtr(blockID,facezData,FACEZ)
endif
#endif
deallocate(xCoord)
deallocate(xCoordL)
deallocate(xCoordR)
deallocate(yCoord)
deallocate(yCoordL)
deallocate(yCoordR)
deallocate(zCoord)
deallocate(zCoordL)
deallocate(zCoordR)
#ifndef FIXEDBLOCKSIZE
deallocate(Az)
deallocate(Ax)
deallocate(Ay)
#endif
end subroutine Simulation_initBlock
############################################################
############################################################
# Simulation_init.F90
############################################################
############################################################
subroutine Simulation_init(myPE)
use Simulation_data
use RuntimeParameters_interface, ONLY : RuntimeParameters_get
use Logfile_interface, ONLY : Logfile_stamp
implicit none
#include "constants.h"
#include "Flash.h"
integer, intent(in) :: myPE
call RuntimeParameters_get('gamma', sim_gamma)
call RuntimeParameters_get('xmin', sim_xMin)
call RuntimeParameters_get('ymin', sim_yMin)
call RuntimeParameters_get('zmin', sim_zMin)
call RuntimeParameters_get('xmax', sim_xMax)
call RuntimeParameters_get('ymax', sim_yMax)
call RuntimeParameters_get('zmax', sim_zMax)
call RuntimeParameters_get('xCtr', sim_xCtr)
call RuntimeParameters_get('yCtr', sim_yCtr)
call RuntimeParameters_get('zCtr', sim_zCtr)
call RuntimeParameters_get('Az_initial', sim_Az_initial)
call RuntimeParameters_get('R_fieldLoop', sim_fieldLoopRadius)
call RuntimeParameters_get('killdivb', sim_killdivb)
call RuntimeParameters_get('smallp', sim_smallP)
call RuntimeParameters_get('smallx', sim_smallX)
call RuntimeParameters_get('sim_rhoLeft', sim_rhoLeft)
call RuntimeParameters_get('sim_rhoRight', sim_rhoRight)
call RuntimeParameters_get('sim_pLeft', sim_pLeft)
call RuntimeParameters_get('sim_pRight', sim_pRight)
call RuntimeParameters_get('sim_uLeft', sim_uLeft)
call RuntimeParameters_get('sim_uRight', sim_uRight)
call RuntimeParameters_get('sim_divide', sim_divide)
sim_gCell = .true.
end subroutine Simulation_init
############################################################
############################################################
# Simulation_initSpecies.F90
############################################################
############################################################
subroutine Simulation_initSpecies()
implicit none
#include "Flash.h"
#include "Multispecies.h"
!! Set up two species, with different plasma parameters
! Similar to Deuterium
call Multispecies_setProperty(DEUT_SPEC, A, 2.0)
call Multispecies_setProperty(DEUT_SPEC, Z, 1.0)
call Multispecies_setProperty(DEUT_SPEC, GAMMA, 1.6666667)
! Similar to Tritium
call Multispecies_setProperty(TRIT_SPEC, A, 3.0)
call Multispecies_setProperty(TRIT_SPEC, Z, 1.0)
call Multispecies_setProperty(TRIT_SPEC, GAMMA, 1.6666667)
end subroutine Simulation_initSpecies
More information about the flash-users
mailing list