Read in the user-supplied parameter-file
| Name | Type | Default | Description |
|---|---|---|---|
| base_filename | character(len=std_len) | None | Base file name for simulation output, which will be followed by a number |
| restart_from_file | character(len=std_len) | None | If not 'unavailable', resume from snapshot with this base file name |
| typefilelog | character(len=std_len) | None | Which type of log to write: 'normal', 'special', 'regression_test' |
| firstprocess | logical | None | If true, call initonegrid_usr upon restarting |
| reset_grid | logical | None | If true, rebuild the AMR grid upon restarting |
| snapshotnext | integer | None | |
| convert | logical | None | If true and restart_from_file is given, convert snapshots to other file formats |
| convert_type | character(len=std_len) | None | Which format to use when converting Options are: tecplot, tecplotCC, vtu, vtuCC, vtuB, vtuBCC, tecplotmpi, tecplotCCmpi, vtumpi, vtuCCmpi, vtuBmpi, vtuBCCmpi, pvtumpi, pvtuCCmpi, pvtuBmpi, pvtuBCCmpi, tecline, teclinempi, onegrid |
| saveprim | logical | None | If true, convert from conservative to primitive variables in output |
| usr_filename | character(len=std_len) | None | User parameter file |
| nwauxio | integer | None | Number of auxiliary variables that are only included in output |
| nocartesian | logical | None | IO switches for conversion |
| w_write | logical | None | |
| writelevel | logical | None | |
| writespshift | double precision | None | |
| length_convert_factor | double precision | None | |
| w_convert_factor | double precision | None | Conversion factors the primitive variables |
| time_convert_factor | double precision | None | Conversion factor for time unit |
| level_io | integer | None | |
| level_io_min | integer | None | |
| level_io_max | integer | None | |
| autoconvert | logical | None | If true, already convert to output format during the run |
| slice_type | character(len=std_len) | None | choose data type of slice: vtu, vtuCC, dat, or csv |
| slicenext | integer | None | |
| collapsenext | integer | None | |
| collapse_type | character(len=std_len) | None | |
| type_endian | integer | None |
| Name | Type | Default | Description |
|---|---|---|---|
| tsave | double precision | None | Save output of type N on times tsave(:, N) |
| itsave | integer | None | Save output of type N on iterations itsave(:, N) |
| dtsave | double precision | None | Repeatedly save output of type N when dtsave(N) simulation time has passed |
| ditsave | integer | None | Repeatedly save output of type N when ditsave(N) time steps have passed |
| nslices | integer | None | Number of slices to output |
| slicedir | integer | None | The slice direction for each slice |
| slicecoord | double precision | None | Slice coordinates, see @ref slices.md |
| collapse | logical | None | If collapse(DIM) is true, generate output integrated over DIM |
| collapseLevel | integer | None | The level at which to produce line-integrated / collapsed output |
| time_between_print | double precision | None | to monitor timeintegration loop at given wall-clock time intervals |
| tsave_log | double precision | None | |
| tsave_dat | double precision | None | |
| tsave_slice | double precision | None | |
| tsave_collapsed | double precision | None | |
| tsave_custom | double precision | None | |
| dtsave_log | double precision | None | |
| dtsave_dat | double precision | None | |
| dtsave_slice | double precision | None | |
| dtsave_collapsed | double precision | None | |
| dtsave_custom | double precision | None | |
| ditsave_log | integer | None | |
| ditsave_dat | integer | None | |
| ditsave_slice | integer | None | |
| ditsave_collapsed | integer | None | |
| ditsave_custom | integer | None | |
| tsavestart_log | double precision | None | |
| tsavestart_dat | double precision | None | |
| tsavestart_slice | double precision | None | |
| tsavestart_collapsed | double precision | None | |
| tsavestart_custom | double precision | None | |
| tsavestart | double precision | None | Start of read out (not counting specified read outs) |
| Name | Type | Default | Description |
|---|---|---|---|
| it_init | integer | None | initial iteration count |
| time_init | double precision | None | Start time for the simulation |
| it_max | integer | None | Stop the simulation after this many time steps have been taken |
| time_max | double precision | None | End time for the simulation |
| dtmin | double precision | None | Stop the simulation when the time step becomes smaller than this value |
| reset_it | logical | None | If true, reset iteration count to 0 |
| reset_time | logical | None | If true, reset iteration count and global_time to original values, and start writing snapshots at index 0 |
| wall_time_max | double precision | None | Ending wall time (in hours) for the simulation |
| final_dt_reduction | logical | None | If true, allow final dt reduction for matching time_max on output |
| Name | Type | Default | Description |
|---|---|---|---|
| time_stepper | character(len=std_len) | None | Which time stepper to use |
| time_integrator | character(len=std_len) | None | Which time integrator to use |
| source_split_usr | logical | None | Use split or unsplit way to add user's source terms, default: unsplit |
| typesourcesplit | character(len=std_len) | None | How to apply dimensional splitting to the source terms, see |
| local_timestep | logical | .false. | each cell has its own timestep or not |
| dimsplit | logical | None | |
| typedimsplit | character(len=std_len) | None | |
| flux_scheme | character(len=std_len) | None | Which flux scheme of spatial discretization to use (per grid level) |
| limiter | character(len=std_len) | None | |
| gradient_limiter | character(len=std_len) | None | |
| cada3_radius | double precision | None | radius of the asymptotic region [0.001, 10], larger means more accurate in smooth region but more overshooting at discontinuities |
| loglimit | logical | None | |
| H_correction | logical | .false. | If true, do H-correction to fix the carbuncle problem at grid-aligned shocks |
| typegrad | character(len=std_len) | None | |
| typediv | character(len=std_len) | None | |
| typecurl | character(len=std_len) | None | type of curl operator |
| flatcd | logical | None | |
| flatsh | logical | None | |
| rk2_alfa | double precision | None | coefficient for rk2_alfa |
| imex222_lambda | double precision | None | IMEX-222(lambda) one-parameter family of schemes |
| ssprk_order | integer | None | SSPRK choice of methods (both threestep and fourstep, Shu-Osher 2N* implementation) also fivestep SSPRK54 |
| rk3_switch | integer | None | RK3 Butcher table |
| imex_switch | integer | None | IMEX_232 choice and parameters |
| schmid_rad1 | double precision | None | |
| schmid_rad2 | double precision | None | |
| schmid_rad3 | double precision | None | |
| phys | character(len=std_len) | None | type of physics to build |
| flux_adaptive_diffusion | logical | .false. | Adaptive LLF diffusion (Rempel et al. 2009): phi in [flux_ad_min, 1] per face per variable |
| flux_ad_min | double precision | 0.0d0 | |
| flux_ad_scale | double precision | 1.0d0 |
| Name | Type | Default | Description |
|---|---|---|---|
| nghostcells | integer | 2 | Number of ghost cells surrounding a grid |
| ghost_copy | logical | .false. | whether copy values instead of interpolation in ghost cells of finer blocks |
| internalboundary | logical | None | if there is an internal boundary |
| typeboundary_min1 | character(len=std_len) | None | |
| typeboundary_min2 | character(len=std_len) | None | |
| typeboundary_min3 | character(len=std_len) | None | |
| typeboundary_max1 | character(len=std_len) | None | |
| typeboundary_max2 | character(len=std_len) | None | |
| typeboundary_max3 | character(len=std_len) | None | |
| save_physical_boundary | logical | None | True for save physical boundary cells in dat files |
| specialboundary | logical | .false. | signal to the precompiler that we need special boundaries: |
| Name | Type | Default | Description |
|---|---|---|---|
| refine_max_level | integer | None | Maximal number of AMR levels |
| nbufferx1 | integer | None | Number of cells as buffer zone |
| nbufferx2 | integer | None | Number of cells as buffer zone |
| nbufferx3 | integer | None | Number of cells as buffer zone |
| refine_threshold | double precision | None | Error tolerance for refinement decision |
| derefine_ratio | double precision | None | |
| refine_criterion | integer | None | |
| amr_wavefilter | double precision | None | |
| max_blocks | integer | 4000 | The maximum number of grid blocks in a processor |
| block_nx1 | integer | None | number of cells for each dimension in grid block excluding ghostcells |
| block_nx2 | integer | None | number of cells for each dimension in grid block excluding ghostcells |
| block_nx3 | integer | None | number of cells for each dimension in grid block excluding ghostcells |
| domain_nx1 | integer | None | number of cells for each dimension in level-one mesh |
| domain_nx2 | integer | None | number of cells for each dimension in level-one mesh |
| domain_nx3 | integer | None | number of cells for each dimension in level-one mesh |
| iprob | integer | None | problem switch allowing different setups in same usr_mod.t |
| xprobmin1 | double precision | None | minimum and maximum domain boundaries for each dimension |
| xprobmin2 | double precision | None | minimum and maximum domain boundaries for each dimension |
| xprobmin3 | double precision | None | minimum and maximum domain boundaries for each dimension |
| xprobmax1 | double precision | None | minimum and maximum domain boundaries for each dimension |
| xprobmax2 | double precision | None | minimum and maximum domain boundaries for each dimension |
| xprobmax3 | double precision | None | minimum and maximum domain boundaries for each dimension |
| w_refine_weight | double precision | None | Weights of variables used to calculate error for mesh refinement |
| prolongprimitive | logical | .false. | |
| coarsenprimitive | logical | .false. | |
| typeprolonglimit | character(len=std_len) | None | Limiter used for prolongation to refined grids and ghost cells |
| logflag | logical | None | |
| tfixgrid | double precision | None | Fix the AMR grid after this time |
| itfixgrid | integer | None | Fix the AMR grid after this many time steps |
| ditregrid | integer | None | Reconstruct the AMR grid once every ditregrid iteration(s) |
| refine_usr | logical | .false. | Specify to use user-defined refinement criterion |
| geometry | character(len=std_len) | 'Cartesian' | coordinate system to build for: 'Cartesian', 'spherical', 'cylindrical', or the logarithmic-radius variants 'logSpherical'/'logCylindrical' (which emit GEOM='spherical'/'cylindrical' plus the LOG_RADIUS flag). Read from &meshlist and turned into the GEOM fypp define by the config system; must be consistent with the set_coordinate_system call in mod_usr. |
| log_r0 | double precision | 0.0d0 | Offset of the logarithmic radial map, r0 in s = ln(1 + r/r0). Only meaningful in a LOG_RADIUS build. Zero - the default - selects the plain map s = ln(r), which cannot reach r = 0. A positive value moves the singularity of the map to r = -r0, so the mesh is uniform (with width r0*ds) for r << r0 and logarithmic for r >> r0, and r = 0 becomes an ordinary face of the grid at s = 0. |
| Name | Type | Default | Description |
|---|---|---|---|
| courantpar | double precision | None | The Courant (CFL) number used for the simulation |
| dtpar | double precision | None | If dtpar is positive, it sets the timestep dt, otherwise courantpar is used to limit the time step based on the Courant condition. |
| dtdiffpar | double precision | None | For resistive MHD, the time step is also limited by the diffusion time: \f$ dt < dtdiffpar \times dx^2/eta \f$ |
| typecourant | character(len=std_len) | None | How to compute the CFL-limited time step. Options are 'maxsum': max(sum(c/dx)); 'summax': sum(max(c/dx)) and 'minimum: max(c/dx), where the summations loop over the grid dimensions and c is the velocity. The default 'maxsum' is the conventiontal way of computing CFL-limited time steps. |
| slowsteps | integer | None | If > 1, then in the first slowsteps-1 time steps dt is reduced by a factor \f$ 1 - (1- step/slowsteps)^2 \f$ |
| Name | Type | Default | Description |
|---|---|---|---|
| filename_euv | character(len=std_len) | None | Base file name for synthetic EUV emission output |
| wavelength | integer | None | wavelength for output |
| filename_sxr | character(len=std_len) | None | Base file name for synthetic SXR emission output |
| emin_sxr | integer | None | |
| emax_sxr | integer | None | |
| LOS_theta | double precision | None | direction of the line of sight (LOS) |
| LOS_phi | double precision | None | direction of the line of sight (LOS) |
| image_rotate | double precision | None | rotation of image |
| x_origin | double precision | None | where the is the origin (X=0,Y=0) of image |
| big_image | logical | None | big image |
| spectrum_wl | integer | None | wave length for spectrum |
| location_slit | double precision | None | location of the slit |
| filename_spectrum | character(len=std_len) | None | Base file name for synthetic EUV spectrum output |
| spectrum_window_min | double precision | None | spectral window |
| spectrum_window_max | double precision | None | spectral window |
| instrument_resolution_factor | double precision | None | times for enhancing spatial resolution for EUV image/spectra |
| activate_unit_arcsec | logical | None | use arcsec as length unit of images/spectra |
| filename_whitelight | character(len=std_len) | None | Base file name for synthetic white light |
| whitelight_instrument | character(len=std_len) | None | white light observation instrument |
| R_occultor | double precision | None | the white light emission below it (unit=Rsun) is not visible |
| R_opt_thick | double precision | None | for spherical coordinate, region below it (unit=Rsun) is treated as not transparent |
| dat_resolution | logical | None | resolution of the images |
| direction_slit | integer | None | direction of the slit (for dat resolution only) |
subroutine read_par_files() use mod_global_parameters use mod_physics, only: phys_energy, physics_type, phys_wider_stencil use mod_limiter use mod_slice use mod_geometry use mod_source use mod_input_output_helper, only: get_names_from_string logical :: fileopen, file_exists integer :: i, j, k, ifile, io_state integer :: iB, isave, iw, level, idim, islice, iside integer :: nx_vec(3), block_nx_vec(3) integer :: my_unit, iostate integer :: ilev double precision :: dx_vec(3) character :: c_ndim character(len=80) :: fmt_string character(len=std_len) :: err_msg, restart_from_file_arg character(len=std_len) :: basename_full, basename_prev, dummy_file character(len=std_len), dimension(:), allocatable :: typeboundary_min1,& typeboundary_min2,typeboundary_min3, typeboundary_max1,& typeboundary_max2,typeboundary_max3 character(len=std_len), allocatable :: limiter(:) character(len=std_len), allocatable :: gradient_limiter(:) !> How to apply dimensional splitting to the source terms, see !> @ref discretization.md character(len=std_len) :: typesourcesplit !> Which flux scheme of spatial discretization to use (per grid level) character(len=std_len), allocatable :: flux_scheme(:) !> Which time stepper to use character(len=std_len) :: time_stepper !> Which time integrator to use character(len=std_len) :: time_integrator !> type of curl operator character(len=std_len) :: typecurl !> Limiter used for prolongation to refined grids and ghost cells character(len=std_len) :: typeprolonglimit !> How to compute the CFL-limited time step. !> Options are 'maxsum': max(sum(c/dx)); 'summax': sum(max(c/dx)) and !> 'minimum: max(c/dx), where the summations loop over the grid dimensions and !> c is the velocity. The default 'maxsum' is the conventiontal way of !> computing CFL-limited time steps. character(len=std_len) :: typecourant double precision, dimension(nsavehi) :: tsave_log, tsave_dat, tsave_slice,& tsave_collapsed, tsave_custom double precision :: dtsave_log, dtsave_dat, dtsave_slice, dtsave_collapsed,& dtsave_custom integer :: ditsave_log, ditsave_dat, ditsave_slice, ditsave_collapsed,& ditsave_custom double precision :: tsavestart_log, tsavestart_dat, tsavestart_slice,& tsavestart_collapsed, tsavestart_custom integer :: ipower double precision :: sizeuniformpart1,sizeuniformpart2,sizeuniformpart3 double precision :: im_delta,im_nu,rka54,rka51,rkb54,rka55 namelist /filelist/ base_filename,restart_from_file, typefilelog,& firstprocess,reset_grid,snapshotnext, convert,convert_type,saveprim,& usr_filename,nwauxio,nocartesian, w_write,writelevel,writespshift,& length_convert_factor, w_convert_factor, time_convert_factor,level_io,& level_io_min, level_io_max, autoconvert,slice_type,slicenext,& collapsenext,collapse_type, type_endian namelist /savelist/ tsave,itsave,dtsave,ditsave,nslices,slicedir,& slicecoord,collapse,collapseLevel, time_between_print,tsave_log,& tsave_dat, tsave_slice, tsave_collapsed, tsave_custom, dtsave_log,& dtsave_dat, dtsave_slice, dtsave_collapsed, dtsave_custom, ditsave_log,& ditsave_dat, ditsave_slice, ditsave_collapsed, ditsave_custom,& tsavestart_log, tsavestart_dat, tsavestart_slice, tsavestart_collapsed,& tsavestart_custom, tsavestart namelist /stoplist/ it_init,time_init,it_max,time_max,dtmin,reset_it,& reset_time,wall_time_max,final_dt_reduction namelist /methodlist/ time_stepper,time_integrator, source_split_usr,& typesourcesplit,local_timestep,dimsplit,typedimsplit,flux_scheme,& limiter,gradient_limiter,cada3_radius,loglimit,& H_correction, typegrad,& typediv,typecurl,flatcd,flatsh,& rk2_alfa,imex222_lambda,ssprk_order,rk3_switch,imex_switch,& schmid_rad1,schmid_rad2,& schmid_rad3, phys,& flux_adaptive_diffusion, flux_ad_min, flux_ad_scale namelist /boundlist/ nghostcells,ghost_copy,internalboundary,& typeboundary_min1,typeboundary_min2,typeboundary_min3,& typeboundary_max1,typeboundary_max2,typeboundary_max3,& save_physical_boundary, specialboundary namelist /meshlist/ refine_max_level,nbufferx1,nbufferx2,nbufferx3,& refine_threshold,derefine_ratio, refine_criterion,& amr_wavefilter, max_blocks,& block_nx1,block_nx2,block_nx3,domain_nx1,& domain_nx2,domain_nx3,iprob,xprobmin1,xprobmin2,xprobmin3,xprobmax1,& xprobmax2,xprobmax3, w_refine_weight, prolongprimitive,coarsenprimitive,& typeprolonglimit, logflag,tfixgrid,itfixgrid,ditregrid, refine_usr,& geometry, log_r0 namelist /paramlist/ courantpar, dtpar, dtdiffpar, typecourant, slowsteps namelist /emissionlist/ filename_euv,wavelength,filename_sxr,emin_sxr,& emax_sxr,LOS_theta,LOS_phi,image_rotate,x_origin,big_image,spectrum_wl,& location_slit,filename_spectrum,spectrum_window_min,spectrum_window_max,& instrument_resolution_factor,activate_unit_arcsec,filename_whitelight,& whitelight_instrument,R_occultor,R_opt_thick,dat_resolution,& direction_slit ! allocate cell size of all levels allocate(dx(ndim,nlevelshi)) allocate(dg1(nlevelshi)) allocate(dg2(nlevelshi)) allocate(dg3(nlevelshi)) allocate(ng1(nlevelshi)) allocate(ng2(nlevelshi)) allocate(ng3(nlevelshi)) ! default block size excluding ghost cells block_nx1 = 16 block_nx2 = 16 block_nx3 = 16 ! default resolution of level-1 mesh (full domain) domain_nx1 = 32 domain_nx2 = 32 domain_nx3 = 32 !! default number of ghost-cell layers at each boundary of a block ! this is now done when the variable is defined in mod_global_parameters ! the physics modules might set this variable in their init subroutine called earlier nghostcells = 2 ! Allocate boundary conditions arrays in new and old style allocate(typeboundary_min1(nwfluxbc)) allocate(typeboundary_max1(nwfluxbc)) typeboundary_min1 = undefined typeboundary_max1 = undefined allocate(typeboundary_min2(nwfluxbc)) allocate(typeboundary_max2(nwfluxbc)) typeboundary_min2 = undefined typeboundary_max2 = undefined allocate(typeboundary_min3(nwfluxbc)) allocate(typeboundary_max3(nwfluxbc)) typeboundary_min3 = undefined typeboundary_max3 = undefined allocate(typeboundary(nwflux+nwaux,2*ndim)) typeboundary=0 ! not save physical boundary in dat files by default save_physical_boundary = .false. internalboundary = .false. ! defaults for specific options typegrad = 'central' typediv = 'central' typecurl = 'central' ! defaults for convert behavior ! store the -if value from argument in command line restart_from_file_arg = restart_from_file nwauxio = 0 nocartesian = .false. saveprim = .false. autoconvert = .false. convert_type = 'vtuBCCmpi' slice_type = 'vtuCC' collapse_type = 'vti' allocate(w_write(nw)) w_write(1:nw) = .true. allocate(writelevel(nlevelshi)) writelevel(1:nlevelshi) = .true. writespshift(1:ndim,1:2) = zero level_io = -1 level_io_min = 1 level_io_max = nlevelshi ! endianness: littleendian (default) is 1, bigendian otherwise type_endian = 1 ! normalization of primitive variables: only for output ! note that length_convert_factor is for length ! this scaling is optional, and must be set consistently if used allocate(w_convert_factor(nw)) w_convert_factor(:) = 1.0d0 time_convert_factor = 1.0d0 length_convert_factor = 1.0d0 ! AMR related defaults refine_max_level = 1 nbufferx1 = 0 nbufferx2 = 0 nbufferx3 = 0 allocate(refine_threshold(nlevelshi)) refine_threshold(1:nlevelshi) = 0.1d0 allocate(derefine_ratio(nlevelshi)) derefine_ratio(1:nlevelshi) = 1.0d0/8.0d0 typeprolonglimit = 'default' refine_criterion = 3 allocate(w_refine_weight(nw+1)) w_refine_weight = 0.d0 allocate(logflag(nw+1)) logflag = .false. allocate(amr_wavefilter(nlevelshi)) amr_wavefilter(1:nlevelshi) = 1.0d-2 tfixgrid = bigdouble itfixgrid = biginteger ditregrid = 1 ! Offset of the logarithmic radial map; 0 selects the plain s = ln(r) log_r0 = 0.0d0 ! IO defaults it_init = 0 it_max = biginteger time_init = 0.d0 time_max = bigdouble wall_time_max = bigdouble if(local_timestep) then final_dt_reduction=.false. else final_dt_reduction=.true. endif final_dt_exit=.false. dtmin = 1.0d-10 nslices = 0 collapse = .false. collapseLevel = 1 time_between_print = 30.0d0 ! Print status every 30 seconds do ifile=1,nfile do isave=1,nsavehi tsave(isave,ifile) = bigdouble !global_time of saves into the output files itsave(isave,ifile) = biginteger ! it of saves into the output files end do dtsave(ifile) = bigdouble ! time between saves ditsave(ifile) = biginteger ! timesteps between saves isavet(ifile) = 1 ! index for saves by global_time isaveit(ifile) = 1 ! index for saves by it tsavestart(ifile) = 0.0d0 end do tsave_log = bigdouble tsave_dat = bigdouble tsave_slice = bigdouble tsave_collapsed = bigdouble tsave_custom = bigdouble dtsave_log = bigdouble dtsave_dat = bigdouble dtsave_slice = bigdouble dtsave_collapsed = bigdouble dtsave_custom = bigdouble ditsave_log = biginteger ditsave_dat = biginteger ditsave_slice = biginteger ditsave_collapsed = biginteger ditsave_custom = biginteger tsavestart_log = bigdouble tsavestart_dat = bigdouble tsavestart_slice = bigdouble tsavestart_collapsed = bigdouble tsavestart_custom = bigdouble typefilelog = 'default' ! defaults for input reset_time = .false. reset_it = .false. firstprocess = .false. reset_grid = .false. base_filename = 'data' usr_filename = '' ! Defaults for discretization methods slowsteps = -1 courantpar = 0.8d0 typecourant = 'maxsum' dimsplit = .false. typedimsplit = 'default' if(physics_type=='mhd') then cada3_radius = 0.5d0 else cada3_radius = 0.1d0 end if schmid_rad1 = 1.d0 schmid_rad2 = 1.d0 schmid_rad3 = 1.d0 source_split_usr= .false. time_stepper = 'twostep' time_integrator = 'default' ! default PC or explicit midpoint, hence alfa=0.5 rk2_alfa = half ! default IMEX-RK22Ln hence lambda = 1 - 1/sqrt(2) imex222_lambda = 1.0d0 - 1.0d0 / dsqrt(2.0d0) ! default SSPRK(3,3) or Gottlieb-Shu 1998 for threestep ! default SSPRK(4,3) or Spireti-Ruuth for fourstep ! default SSPRK(5,4) using Gottlieb coeffs ssprk_order = 3 ! default RK3 butcher table: Heun 3rd order rk3_switch = 3 ! default IMEX threestep is IMEX_ARK(232) imex_switch = 1 ! Defaults for synthesing emission LOS_theta = 0.d0 LOS_phi = 0.d0 image_rotate = 0.d0 x_origin = 0.d0 big_image = .false. location_slit = 0.d0 direction_slit = -1 instrument_resolution_factor=1.d0 activate_unit_arcsec=.true. whitelight_instrument='LASCO/C2' R_occultor=-1.d0 R_opt_thick=1.d0 dat_resolution=.false. allocate(flux_scheme(nlevelshi), flux_method(nlevelshi)) allocate(limiter(nlevelshi),gradient_limiter(nlevelshi)) do level=1,nlevelshi flux_scheme(level) = 'tvdlf' limiter(level) = 'minmod' gradient_limiter(level) = 'minmod' end do flatcd = .false. flatsh = .false. typesourcesplit = 'sfs' allocate(loglimit(nw)) loglimit(1:nw) = .false. dtdiffpar = 0.5d0 dtpar = -1.d0 ! problem setup defaults iprob = 1 ! end defaults ! Initialize Kronecker delta, and Levi-Civita tensor do i=1,3 do j=1,3 if(i==j)then kr(i,j)=1 else kr(i,j)=0 endif do k=1,3 if(i==j.or.j==k.or.k==i)then lvc(i,j,k)=0 else if(i+1==j.or.i-2==j)then lvc(i,j,k)=1 else lvc(i,j,k)=-1 endif enddo enddo enddo ! These are used to construct file and log names from multiple par files basename_full = '' basename_prev = '' ! Cray compiler bug: setting values after reading the namelists fails ! when reading files in a loop. ! To avoid this, disallow reading more than one file for now #ifdef _CRAYFTN if (size(par_files) > 1) then write(err_msg, *) "Reading more than one par file is disallowed " // & "when compiled with Cray to avoid a compiler bug" call mpistop(trim(err_msg)) end if if (size(par_files) > 0) then i = 1 #else do i = 1, size(par_files) #endif if (mype == 0) print *, "Reading " // trim(par_files(i)) ! Check whether the file exists inquire(file=trim(par_files(i)), exist=file_exists) if (.not. file_exists) then write(err_msg, *) "The parameter file " // trim(par_files(i)) // & " does not exist" call mpistop(trim(err_msg)) end if open(unitpar, file=trim(par_files(i)), status='old') ! Try to read in the namelists. They can be absent or in a different ! order, since we rewind before each read. rewind(unitpar) read(unitpar, filelist, end=101) 101 rewind(unitpar) read(unitpar, savelist, end=102) 102 rewind(unitpar) read(unitpar, stoplist, end=103) 103 rewind(unitpar) read(unitpar, methodlist, end=104) 104 rewind(unitpar) read(unitpar, boundlist, end=105) 105 rewind(unitpar) read(unitpar, meshlist, end=106) 106 rewind(unitpar) read(unitpar, paramlist, end=107) 107 rewind(unitpar) read(unitpar, emissionlist, end=108) 108 close(unitpar) ! Append the log and file names given in the par files if (base_filename /= basename_prev) basename_full = trim(basename_full) & // trim(base_filename) basename_prev = base_filename #ifdef _CRAYFTN end if #else end do #endif base_filename = basename_full ! Check whether output directory is writable if(mype==0) then dummy_file = trim(base_filename)//"DUMMY" open(newunit=my_unit, file=trim(dummy_file), iostat=iostate) if (iostate /= 0) then call mpistop("Can't write to output directory (" // & trim(base_filename) // ")") else close(my_unit, status='delete') end if end if if(source_split_usr) any_source_split=.true. ! restart filename from command line overwrites the one in par file if(restart_from_file_arg /= undefined) & restart_from_file=restart_from_file_arg ! Root process will search snapshot if (mype == 0) then if(restart_from_file == undefined) then ! search file from highest index file_exists=.false. do index_latest_data = 9999, 0, -1 if(snapshot_exists(index_latest_data)) then file_exists=.true. exit end if end do if(.not.file_exists) index_latest_data=-1 else ! get index of the given data restarted from index_latest_data=get_snapshot_index(trim(restart_from_file)) end if end if call MPI_BCAST(index_latest_data, 1, MPI_INTEGER, 0, icomm, ierrmpi) if (resume_previous_run) then if (index_latest_data == -1) then if(mype==0) write(*,*)& "No snapshots found to resume from, start a new run..." else ! Set file name to restart from write(restart_from_file, "(a,i4.4,a)") trim(base_filename),& index_latest_data, ".dat" end if end if if (restart_from_file == undefined) then snapshotnext = 0 slicenext = 0 collapsenext = 0 if (firstprocess) call mpistop& ("Please restart from a snapshot when firstprocess=T") if (convert) call mpistop('Change convert to .false. for a new run!') end if if(convert) autoconvert=.false. where (tsave_log < bigdouble) tsave(:, 1) = tsave_log where (tsave_dat < bigdouble) tsave(:, 2) = tsave_dat where (tsave_slice < bigdouble) tsave(:, 3) = tsave_slice where (tsave_collapsed < bigdouble) tsave(:, 4) = tsave_collapsed where (tsave_custom < bigdouble) tsave(:, 5) = tsave_custom if (dtsave_log < bigdouble) dtsave(1) = dtsave_log if (dtsave_dat < bigdouble) dtsave(2) = dtsave_dat if (dtsave_slice < bigdouble) dtsave(3) = dtsave_slice if (dtsave_collapsed < bigdouble) dtsave(4) = dtsave_collapsed if (dtsave_custom < bigdouble) dtsave(5) = dtsave_custom if (tsavestart_log < bigdouble) tsavestart(1) = tsavestart_log if (tsavestart_dat < bigdouble) tsavestart(2) = tsavestart_dat if (tsavestart_slice < bigdouble) tsavestart(3) = tsavestart_slice if (tsavestart_collapsed < bigdouble) tsavestart(4) = tsavestart_collapsed if (tsavestart_custom < bigdouble) tsavestart(5) = tsavestart_custom if (ditsave_log < bigdouble) ditsave(1) = ditsave_log if (ditsave_dat < bigdouble) ditsave(2) = ditsave_dat if (ditsave_slice < bigdouble) ditsave(3) = ditsave_slice if (ditsave_collapsed < bigdouble) ditsave(4) = ditsave_collapsed if (ditsave_custom < bigdouble) ditsave(5) = ditsave_custom ! convert hours to seconds for ending wall time if (wall_time_max < bigdouble) wall_time_max=wall_time_max*3600.d0 if (mype == 0) then write(unitterm, *) '' write(unitterm, *)& 'Output type | tsavestart | dtsave | ditsave | itsave(1) | tsave(1)' write(fmt_string, *) '(A12," | ",E9.3E2," | ",E9.3E2," | ",I6," | "'//& ',I6, " | ",E9.3E2)' end if do ifile = 1, nfile if (mype == 0) write(unitterm, fmt_string) trim(output_names(ifile)),& tsavestart(ifile), dtsave(ifile), ditsave(ifile), itsave(1, ifile),& tsave(1, ifile) end do if (mype == 0) write(unitterm, *) '' do islice=1,nslices if(slicedir(islice) > ndim) write(uniterr,& *)'Warning in read_par_files: ', 'Slice ', islice,' direction',& slicedir(islice),'larger than ndim=',ndim if(slicedir(islice) < 1) write(uniterr,*)'Warning in read_par_files: ',& 'Slice ', islice,' direction',slicedir(islice),& 'too small, should be [',1,ndim,']' end do if(it_max==biginteger .and. time_max==bigdouble.and.mype==0) write(uniterr,& *) 'Warning in read_par_files: it_max or time_max not given!' select case (typecourant) case ('maxsum') type_courant=type_maxsum case ('summax') type_courant=type_summax case ('minimum') type_courant=type_minimum case default write(unitterm,*)'Unknown typecourant=',typecourant call mpistop("Error from read_par_files: no such typecourant!") end select ! ========================================================================= ! Flux scheme parsing and physics compatibility validation ! ========================================================================= do level=1,nlevelshi select case (flux_scheme(level)) case ('hll') flux_method(level)=fs_hll if (physics_type=='ffhd') then call mpistop(trim(flux_scheme(level))//" flux scheme is not compatible with "//trim(physics_type)) end if case ('hllc') flux_method(level)=fs_hllc if (physics_type=='mhd' .or. physics_type=='ffhd') then call mpistop(trim(flux_scheme(level))//" flux scheme is not compatible with "//trim(physics_type)) end if case ('tvdlf') flux_method(level)=fs_tvdlf case default call mpistop("unknown or bad flux scheme: "//trim(flux_scheme(level))) end select end do ! ========================================================================= ! Limiter and spatial operator parsing ! ========================================================================= if(any(limiter=='schmid1')) need_global_a2max=.true. ! initialize type_curl select case (typecurl) case ("central") type_curl=central case ("Gaussbased") type_curl=Gaussbased case ("Stokesbased") type_curl=Stokesbased case default write(unitterm,*) "typecurl=",typecurl call mpistop("unkown type of curl operator in read_par_files") end select ! initialize types of time stepper and time integrator select case (time_stepper) case ("onestep") t_stepper=onestep nstep=1 if (time_integrator=='default') then time_integrator="Forward_Euler" end if select case (time_integrator) case ("Forward_Euler") t_integrator=Forward_Euler case ("IMEX_Euler") t_integrator=IMEX_Euler case ("IMEX_SP") t_integrator=IMEX_SP case default write(unitterm,*) "time_integrator=",time_integrator,"time_stepper=",& time_stepper call mpistop("unkown onestep time_integrator in read_par_files") end select use_imex_scheme=(t_integrator==IMEX_Euler.or.t_integrator==IMEX_SP) case ("twostep") t_stepper=twostep nstep=2 if (time_integrator=='default') then time_integrator="Predictor_Corrector" endif select case (time_integrator) case ("Predictor_Corrector") t_integrator=Predictor_Corrector case ("RK2_alfa") t_integrator=RK2_alf case ("ssprk2") t_integrator=ssprk2 case ("IMEX_Midpoint") t_integrator=IMEX_Midpoint case ("IMEX_Trapezoidal") t_integrator=IMEX_Trapezoidal case ("IMEX_222") t_integrator=IMEX_222 case default write(unitterm,*) "time_integrator=",time_integrator,"time_stepper=",& time_stepper call mpistop("unkown twostep time_integrator in read_par_files") end select use_imex_scheme=(t_integrator==IMEX_Midpoint.or.t_integrator==& IMEX_Trapezoidal.or.t_integrator==IMEX_222) if (t_integrator==RK2_alf) then if(rk2_alfa<smalldouble.or.rk2_alfa>one)call & mpistop("set rk2_alfa within [0,1]") rk_a21=rk2_alfa rk_b2=1.0d0/(2.0d0*rk2_alfa) rk_b1=1.0d0-rk_b2 endif case ("threestep") t_stepper=threestep nstep=3 if (time_integrator=='default') then time_integrator='ssprk3' endif select case (time_integrator) case ("ssprk3") t_integrator=ssprk3 case ("RK3_BT") t_integrator=RK3_BT case ("IMEX_ARS3") t_integrator=IMEX_ARS3 case ("IMEX_232") t_integrator=IMEX_232 case ("IMEX_CB3a") t_integrator=IMEX_CB3a case default write(unitterm,*) "time_integrator=",time_integrator,"time_stepper=",& time_stepper call mpistop("unkown threestep time_integrator in read_par_files") end select if(t_integrator==RK3_BT) then select case(rk3_switch) case(1) ! we code up Ralston 3rd order here rk3_a21=1.0d0/2.0d0 rk3_a31=0.0d0 rk3_a32=3.0d0/4.0d0 rk3_b1=2.0d0/9.0d0 rk3_b2=1.0d0/3.0d0 case(2) ! we code up RK-Wray 3rd order here rk3_a21=8.0d0/15.0d0 rk3_a31=1.0d0/4.0d0 rk3_a32=5.0d0/12.0d0 rk3_b1=1.0d0/4.0d0 rk3_b2=0.0d0 case(3) ! we code up Heun 3rd order here rk3_a21=1.0d0/3.0d0 rk3_a31=0.0d0 rk3_a32=2.0d0/3.0d0 rk3_b1=1.0d0/4.0d0 rk3_b2=0.0d0 case(4) ! we code up Nystrom 3rd order here rk3_a21=2.0d0/3.0d0 rk3_a31=0.0d0 rk3_a32=2.0d0/3.0d0 rk3_b1=1.0d0/4.0d0 rk3_b2=3.0d0/8.0d0 case default call mpistop("Unknown rk3_switch") end select ! the rest is fixed from above rk3_b3=1.0d0-rk3_b1-rk3_b2 rk3_c2=rk3_a21 rk3_c3=rk3_a31+rk3_a32 endif if(t_integrator==ssprk3) then select case(ssprk_order) case(3) ! this is SSPRK(3,3) Gottlieb-Shu rk_beta11=1.0d0 rk_beta22=1.0d0/4.0d0 rk_beta33=2.0d0/3.0d0 rk_alfa21=3.0d0/4.0d0 rk_alfa31=1.0d0/3.0d0 rk_c2=1.0d0 rk_c3=1.0d0/2.0d0 case(2) ! this is SSP(3,2) rk_beta11=1.0d0/2.0d0 rk_beta22=1.0d0/2.0d0 rk_beta33=1.0d0/3.0d0 rk_alfa21=0.0d0 rk_alfa31=1.0d0/3.0d0 rk_c2=1.0d0/2.0d0 rk_c3=1.0d0 case default call mpistop("Unknown ssprk3_order") end select rk_alfa22=1.0d0-rk_alfa21 rk_alfa33=1.0d0-rk_alfa31 endif if(t_integrator==IMEX_ARS3) then ars_gamma=(3.0d0+dsqrt(3.0d0))/6.0d0 endif if(t_integrator==IMEX_232) then select case(imex_switch) case(1) ! this is IMEX_ARK(232) im_delta=1.0d0-1.0d0/dsqrt(2.0d0) im_nu=(3.0d0+2.0d0*dsqrt(2.0d0))/6.0d0 imex_a21=2.0d0*im_delta imex_a31=1.0d0-im_nu imex_a32=im_nu imex_b1=1.0d0/(2.0d0*dsqrt(2.0d0)) imex_b2=1.0d0/(2.0d0*dsqrt(2.0d0)) imex_ha21=im_delta imex_ha22=im_delta case(2) ! this is IMEX_SSP(232) ! doi 10.1002/2017MS001065 Rokhzadi et al imex_a21=0.711664700366941d0 imex_a31=0.077338168947683d0 imex_a32=0.917273367886007d0 imex_b1=0.398930808264688d0 imex_b2=0.345755244189623d0 imex_ha21=0.353842865099275d0 imex_ha22=0.353842865099275d0 case default call mpistop("Unknown imex_siwtch") end select imex_c2=imex_a21 imex_c3=imex_a31+imex_a32 imex_b3=1.0d0-imex_b1-imex_b2 endif if(t_integrator==IMEX_CB3a) then imex_c2 = 0.8925502329346865 imex_a22 = imex_c2 imex_ha21 = imex_c2 imex_c3 = imex_c2 / (6.0d0*imex_c2**2 - 3.0d0*imex_c2 + 1.0d0) imex_ha32 = imex_c3 imex_b2 = (3.0d0*imex_c2 - 1.0d0) / (6.0d0*imex_c2**2) imex_b3 = (6.0d0*imex_c2**2 - 3.0d0*imex_c2 + 1.0d0) / & (6.0d0*imex_c2**2) imex_a33 = (1.0d0/6.0d0 - imex_b2*imex_c2**2 - & imex_b3*imex_c2*imex_c3) / (imex_b3*(imex_c3-imex_c2)) imex_a32 = imex_c3 - imex_a33 ! if (mype == 0) then ! write(*,*) "=================================" ! write(*,*) "Asserting the order conditions..." ! ! First order convergence: OK ! ! Second order convergence ! write(*,*) -1.0d0/2.0d0 + imex_b2*imex_c2 + imex_b3*imex_c3 ! ! Third order convergence ! write(*,*) -1.0d0/3.0d0 + imex_b2*imex_c2**2 + imex_b3*imex_c3**2 ! write(*,*) -1.0d0/6.0d0 + imex_b3*imex_ha32*imex_c2 ! write(*,*) -1.0d0/6.0d0 + imex_b2*imex_a22*imex_c2 + imex_b3*imex_a32*imex_c2 + imex_b3*imex_a33*imex_c3 ! write(*,*) "=================================" ! end if end if use_imex_scheme=(t_integrator==IMEX_ARS3.or.t_integrator==& IMEX_232.or.t_integrator==IMEX_CB3a) case ("fourstep") t_stepper=fourstep nstep=4 if (time_integrator=='default') then time_integrator="ssprk4" endif select case (time_integrator) case ("ssprk4") t_integrator=ssprk4 case ("rk4") t_integrator=rk4 case default write(unitterm,*) "time_integrator=",time_integrator,"time_stepper=",& time_stepper call mpistop("unkown fourstep time_integrator in read_par_files") end select if(t_integrator==ssprk4) then select case(ssprk_order) case(3) ! this is SSPRK(4,3) Spireti-Ruuth rk_beta11=1.0d0/2.0d0 rk_beta22=1.0d0/2.0d0 rk_beta33=1.0d0/6.0d0 rk_beta44=1.0d0/2.0d0 rk_alfa21=0.0d0 rk_alfa31=2.0d0/3.0d0 rk_alfa41=0.0d0 rk_c2=1.0d0/2.0d0 rk_c3=1.0d0 rk_c4=1.0d0/2.0d0 case(2) ! this is SSP(4,2) rk_beta11=1.0d0/3.0d0 rk_beta22=1.0d0/3.0d0 rk_beta33=1.0d0/3.0d0 rk_beta44=1.0d0/4.0d0 rk_alfa21=0.0d0 rk_alfa31=0.0d0 rk_alfa41=1.0d0/4.0d0 rk_c2=1.0d0/3.0d0 rk_c3=2.0d0/3.0d0 rk_c4=1.0d0 case default call mpistop("Unknown ssprk_order") end select rk_alfa22=1.0d0-rk_alfa21 rk_alfa33=1.0d0-rk_alfa31 rk_alfa44=1.0d0-rk_alfa41 endif case ("fivestep") t_stepper=fivestep nstep=5 if (time_integrator=='default') then time_integrator="ssprk5" end if select case (time_integrator) case ("ssprk5") t_integrator=ssprk5 case default write(unitterm,*) "time_integrator=",time_integrator,"time_stepper=",& time_stepper call mpistop("unkown fivestep time_integrator in read_par_files") end select if(t_integrator==ssprk5) then select case(ssprk_order) ! we use ssprk_order to intercompare the different coefficient choices case(3) ! From Gottlieb 2005 rk_beta11=0.391752226571890d0 rk_beta22=0.368410593050371d0 rk_beta33=0.251891774271694d0 rk_beta44=0.544974750228521d0 rk_beta54=0.063692468666290d0 rk_beta55=0.226007483236906d0 rk_alfa21=0.444370493651235d0 rk_alfa31=0.620101851488403d0 rk_alfa41=0.178079954393132d0 rk_alfa53=0.517231671970585d0 rk_alfa54=0.096059710526147d0 rk_alfa22=0.555629506348765d0 rk_alfa33=0.379898148511597d0 rk_alfa44=0.821920045606868d0 rk_alfa55=0.386708617503269d0 rk_alfa22=1.0d0-rk_alfa21 rk_alfa33=1.0d0-rk_alfa31 rk_alfa44=1.0d0-rk_alfa41 rk_alfa55=1.0d0-rk_alfa53-rk_alfa54 rk_c2=rk_beta11 rk_c3=rk_alfa22*rk_c2+rk_beta22 rk_c4=rk_alfa33*rk_c3+rk_beta33 rk_c5=rk_alfa44*rk_c4+rk_beta44 case(2) ! From Spireti-Ruuth rk_beta11=0.39175222700392d0 rk_beta22=0.36841059262959d0 rk_beta33=0.25189177424738d0 rk_beta44=0.54497475021237d0 !rk_beta54=0.06369246925946d0 !rk_beta55=0.22600748319395d0 rk_alfa21=0.44437049406734d0 rk_alfa31=0.62010185138540d0 rk_alfa41=0.17807995410773d0 rk_alfa53=0.51723167208978d0 !rk_alfa54=0.09605971145044d0 rk_alfa22=1.0d0-rk_alfa21 rk_alfa33=1.0d0-rk_alfa31 rk_alfa44=1.0d0-rk_alfa41 rka51=0.00683325884039d0 rka54=0.12759831133288d0 rkb54=0.08460416338212d0 rk_beta54=rkb54-rk_beta44*rka51/rk_alfa41 rk_alfa54=rka54-rk_alfa44*rka51/rk_alfa41 rk_alfa55=1.0d0-rk_alfa53-rk_alfa54 rk_c2=rk_beta11 rk_c3=rk_alfa22*rk_c2+rk_beta22 rk_c4=rk_alfa33*rk_c3+rk_beta33 rk_c5=rk_alfa44*rk_c4+rk_beta44 rk_beta55=1.0d0-rk_beta54-rk_alfa53*rk_c3-rk_alfa54*rk_c4-& rk_alfa55*rk_c5 case default call mpistop("Unknown ssprk_order") end select ! the following combinations must be unity !print *,rk_beta55+rk_beta54+rk_alfa53*rk_c3+rk_alfa54*rk_c4+rk_alfa55*rk_c5 !print *,rk_alfa22+rk_alfa21 !print *,rk_alfa33+rk_alfa31 !print *,rk_alfa44+rk_alfa41 !print *,rk_alfa55+rk_alfa53+rk_alfa54 endif use_imex_scheme=.false. case default call mpistop("Unknown time_stepper in read_par_files") end select ! Harmonize the parameters for dimensional splitting and source splitting if(typedimsplit =='default'.and. dimsplit) typedimsplit='xyyx' if(typedimsplit =='default'.and..not.dimsplit) typedimsplit='unsplit' dimsplit = typedimsplit /='unsplit' ! initialize types of split-source addition select case (typesourcesplit) case ('sfs') sourcesplit=sourcesplit_sfs case ('sf') sourcesplit=sourcesplit_sf case ('ssf') sourcesplit=sourcesplit_ssf case ('ssfss') sourcesplit=sourcesplit_ssfss case default write(unitterm,*)'No such typesourcesplit=',typesourcesplit call mpistop("Error: Unknown typesourcesplit!") end select if(coordinate==-1) then coordinate=Cartesian if(mype==0) then write(*,*) 'Warning: coordinate system is not specified!' write(*,*) 'call set_coordinate_system in usr_init in mod_usr.t' write(*,*) 'Now use Cartesian coordinate' end if end if if(coordinate==Cartesian) then slab=.true. slab_uniform=.true. else slab=.false. slab_uniform=.false. end if if(coordinate==spherical) then if(dimsplit) then if(mype==0)print *& ,'Warning: spherical symmetry needs dimsplit=F, resetting' dimsplit=.false. end if end if if (ndim==1) dimsplit=.false. ! type limiter of prolongation select case(typeprolonglimit) case('unlimit') ! unlimited prolong_limiter=1 case('minmod') prolong_limiter=2 case('woodward') prolong_limiter=3 case('koren') prolong_limiter=4 case default prolong_limiter=0 end select ! Type limiter is of integer type for performance allocate(type_limiter(nlevelshi)) allocate(type_gradient_limiter(nlevelshi)) do level=1,nlevelshi type_limiter(level) = limiter_type(limiter(level)) type_gradient_limiter(level) = limiter_type(gradient_limiter(level)) end do if (any(limiter(1:nlevelshi)== 'ppm').and.(flatsh.and.physics_type==& 'rho')) then call mpistop& (" PPM with flatsh=.true. can not be used with physics_type='rho'!") end if ! Copy boundary conditions to typeboundary, which is used internally do iw=1,nwfluxbc select case(typeboundary_min1(iw)) case("special") typeboundary(iw,2*1-1)=bc_special case("cont") typeboundary(iw,2*1-1)=bc_cont case("symm") typeboundary(iw,2*1-1)=bc_symm case("asymm") typeboundary(iw,2*1-1)=bc_asymm case("periodic") typeboundary(iw,2*1-1)=bc_periodic case("aperiodic") typeboundary(iw,2*1-1)=bc_aperiodic case("noinflow") typeboundary(iw,2*1-1)=bc_noinflow case("pole") typeboundary(iw,2*1-1)=12 case("bc_data") typeboundary(iw,2*1-1)=bc_data case("bc_icarus") typeboundary(iw,2*1-1)=bc_icarus case("character") typeboundary(iw,2*1-1)=bc_character case default write (unitterm,*) "Undefined boundarytype found in read_par_files",& typeboundary_min1(iw),"for variable iw=",iw," and side iB=",2*1-1 end select end do do iw=1,nwfluxbc select case(typeboundary_max1(iw)) case("special") typeboundary(iw,2*1)=bc_special case("cont") typeboundary(iw,2*1)=bc_cont case("symm") typeboundary(iw,2*1)=bc_symm case("asymm") typeboundary(iw,2*1)=bc_asymm case("periodic") typeboundary(iw,2*1)=bc_periodic case("aperiodic") typeboundary(iw,2*1)=bc_aperiodic case("noinflow") typeboundary(iw,2*1)=bc_noinflow case("pole") typeboundary(iw,2*1)=12 case("bc_data") typeboundary(iw,2*1)=bc_data case("bc_icarus") typeboundary(iw,2*1-1)=bc_icarus case("bc_character") typeboundary(iw,2*1)=bc_character case default write (unitterm,*) "Undefined boundarytype found in read_par_files",& typeboundary_max1(iw),"for variable iw=",iw," and side iB=",2*1 end select end do do iw=1,nwfluxbc select case(typeboundary_min2(iw)) case("special") typeboundary(iw,2*2-1)=bc_special case("cont") typeboundary(iw,2*2-1)=bc_cont case("symm") typeboundary(iw,2*2-1)=bc_symm case("asymm") typeboundary(iw,2*2-1)=bc_asymm case("periodic") typeboundary(iw,2*2-1)=bc_periodic case("aperiodic") typeboundary(iw,2*2-1)=bc_aperiodic case("noinflow") typeboundary(iw,2*2-1)=bc_noinflow case("pole") typeboundary(iw,2*2-1)=12 case("bc_data") typeboundary(iw,2*2-1)=bc_data case("bc_icarus") typeboundary(iw,2*2-1)=bc_icarus case("character") typeboundary(iw,2*2-1)=bc_character case default write (unitterm,*) "Undefined boundarytype found in read_par_files",& typeboundary_min2(iw),"for variable iw=",iw," and side iB=",2*2-1 end select end do do iw=1,nwfluxbc select case(typeboundary_max2(iw)) case("special") typeboundary(iw,2*2)=bc_special case("cont") typeboundary(iw,2*2)=bc_cont case("symm") typeboundary(iw,2*2)=bc_symm case("asymm") typeboundary(iw,2*2)=bc_asymm case("periodic") typeboundary(iw,2*2)=bc_periodic case("aperiodic") typeboundary(iw,2*2)=bc_aperiodic case("noinflow") typeboundary(iw,2*2)=bc_noinflow case("pole") typeboundary(iw,2*2)=12 case("bc_data") typeboundary(iw,2*2)=bc_data case("bc_icarus") typeboundary(iw,2*2-1)=bc_icarus case("bc_character") typeboundary(iw,2*2)=bc_character case default write (unitterm,*) "Undefined boundarytype found in read_par_files",& typeboundary_max2(iw),"for variable iw=",iw," and side iB=",2*2 end select end do do iw=1,nwfluxbc select case(typeboundary_min3(iw)) case("special") typeboundary(iw,2*3-1)=bc_special case("cont") typeboundary(iw,2*3-1)=bc_cont case("symm") typeboundary(iw,2*3-1)=bc_symm case("asymm") typeboundary(iw,2*3-1)=bc_asymm case("periodic") typeboundary(iw,2*3-1)=bc_periodic case("aperiodic") typeboundary(iw,2*3-1)=bc_aperiodic case("noinflow") typeboundary(iw,2*3-1)=bc_noinflow case("pole") typeboundary(iw,2*3-1)=12 case("bc_data") typeboundary(iw,2*3-1)=bc_data case("bc_icarus") typeboundary(iw,2*3-1)=bc_icarus case("character") typeboundary(iw,2*3-1)=bc_character case default write (unitterm,*) "Undefined boundarytype found in read_par_files",& typeboundary_min3(iw),"for variable iw=",iw," and side iB=",2*3-1 end select end do do iw=1,nwfluxbc select case(typeboundary_max3(iw)) case("special") typeboundary(iw,2*3)=bc_special case("cont") typeboundary(iw,2*3)=bc_cont case("symm") typeboundary(iw,2*3)=bc_symm case("asymm") typeboundary(iw,2*3)=bc_asymm case("periodic") typeboundary(iw,2*3)=bc_periodic case("aperiodic") typeboundary(iw,2*3)=bc_aperiodic case("noinflow") typeboundary(iw,2*3)=bc_noinflow case("pole") typeboundary(iw,2*3)=12 case("bc_data") typeboundary(iw,2*3)=bc_data case("bc_icarus") typeboundary(iw,2*3-1)=bc_icarus case("bc_character") typeboundary(iw,2*3)=bc_character case default write (unitterm,*) "Undefined boundarytype found in read_par_files",& typeboundary_max3(iw),"for variable iw=",iw," and side iB=",2*3 end select end do ! psi, tracers take the same boundary type as the first variable if (nwfluxbc<nwflux) then do iw = nwfluxbc+1, nwflux typeboundary(iw,:) = typeboundary(1, :) end do end if ! auxiliary variables take the same boundary type as the first variable if (nwaux>0) then do iw = nwflux+1, nwflux+nwaux typeboundary(iw,:) = typeboundary(1, :) end do end if if (any(typeboundary == 0)) then call mpistop("Not all boundary conditions have been defined") end if do idim=1,ndim periodB(idim)=(any(typeboundary(:,2*idim-1:2*idim)==bc_periodic)) aperiodB(idim)=(any(typeboundary(:,2*idim-1:2*idim)==bc_aperiodic)) if (periodB(idim).or.aperiodB(idim)) then do iw=1,nwflux if (typeboundary(iw,2*idim-1) .ne. typeboundary(iw,& 2*idim)) call mpistop("Wrong counterpart in periodic boundary") if (typeboundary(iw,2*idim-1) /= bc_periodic .and. & typeboundary(iw,2*idim-1) /= bc_aperiodic) then call mpistop("Each dimension should either have all "//& "or no variables periodic, some can be aperiodic") end if end do end if end do ! Expand the 'pole' boundary type (parsed as the sentinel 12 above) into ! ordinary per-variable symm/asymm entries. The axis is not a boundary: a ! ghost cell beyond it is the interior cell half a revolution away in phi, ! and getbc fetches it from there. All this table has to say is which ! components are odd under that pi-rotation, and the rule is the same in ! both coordinate systems: the component along the mirrored direction ! flips, the phi component flips, everything else is symmetric. ! ! spherical, axis at theta=0 or pi (mirror in theta): m_theta, m_phi ! cylindrical, axis at r=0 (mirror in r): m_r, m_phi ! ! and likewise for the magnetic field in an mhd build. Note the ! cylindrical radial entry: upstream MPI-AMRVAC marks only m_phi (and ! B_phi) antisymmetric at the cylindrical axis, but the axis maps ! (r, phi) -> (r, phi+pi), under which e_r(phi+pi) = -e_r(phi) exactly as ! e_phi(phi+pi) = -e_phi(phi), so the radial component is odd too. ! tests/hd/cylindrical_pole_uniform_flow is what pins this down. ! ! The indices come from iw_mom/iw_mag rather than from the positional ! rho-mom-[e]-B arithmetic upstream uses, because a HYPERTC build ! registers q_ between e_ and mag and shifts every field index. do idim=1,ndim do iside=1,2 iB=2*idim-2+iside if(.not.any(typeboundary(:,iB)==12)) cycle poleB_requested(iside,idim)=.true. ! a pole is a property of the face, not of one variable if(any(typeboundary(:,iB)/=12)) call mpistop("typeboundary='pole' has & &to be given for every variable on a face") select case(physics_type) case('hd','mhd','srhd') ! the vector components handled below are all there is to say. srhd's ! mom(:) is the spatial four-velocity u^i = lfac*v^i (see con2prim); ! lfac is a scalar that a coordinate rotation leaves invariant, so ! u^i transforms exactly like an ordinary vector under the pole's ! pi-rotation and takes the same asymm pattern as hd's momentum. case('ffhd') ! ffhd's only conserved quantities are scalars - rho, the ! field-aligned momentum mom(1)=m_par, the energy - and a scalar is ! even under the pole's proper pi-rotation (m_par = rho*(v.bhat) is a ! contraction of two vectors), so every getbc-exchanged variable ! stays symm and the momentum lines below are skipped. The frozen ! field b1,b2,b3 is a vector, but it is not exchanged through getbc: ! fill_nwextra_device re-derives it in every ghost cell, pole ! ghosts included, so it needs no entry here either. case default call mpistop('pole treatment is only implemented for phys=hd, & &phys=mhd, phys=srhd and phys=ffhd') end select select case(coordinate) case(cylindrical) if(idim/=1) call mpistop('the cylindrical axis is at r=0, so only & &typeboundary_min1 can be pole') case(spherical) if(idim/=2) call mpistop('the spherical poles are at theta=0 and & &theta=pi, so only typeboundary_min2/max2 can be pole') case default call mpistop('Pole is in cylindrical, polar, spherical coordinates!') end select typeboundary(:,iB)=bc_symm if(physics_type/='ffhd') then ! ffhd's mom has a single element (the scalar m_par); hd/mhd/srhd ! carry a genuine ndir-vector whose mirrored and phi components flip typeboundary(iw_mom(idim),iB)=bc_asymm typeboundary(iw_mom(phi_),iB)=bc_asymm end if if(physics_type=='mhd') then typeboundary(iw_mag(idim),iB)=bc_asymm typeboundary(iw_mag(phi_),iB)=bc_asymm end if end do end do if (any(typeboundary(:,:)==bc_special) .and. .not. specialboundary) then call mpistop('special boundary requested, set specialboundary=.true.') end if if(.not.phys_energy) then flatcd=.false. flatsh=.false. end if if(any(limiter(1:nlevelshi)=='mp5')) then nghostcells=max(nghostcells,3) end if if(any(limiter(1:nlevelshi)=='weno5')) then nghostcells=max(nghostcells,3) end if if(any(limiter(1:nlevelshi)=='weno5nm')) then nghostcells=max(nghostcells,3) end if if(any(limiter(1:nlevelshi)=='wenoz5')) then nghostcells=max(nghostcells,3) end if if(any(limiter(1:nlevelshi)=='wenoz5nm')) then nghostcells=max(nghostcells,3) end if if(any(limiter(1:nlevelshi)=='wenozp5')) then nghostcells=max(nghostcells,3) end if if(any(limiter(1:nlevelshi)=='wenozp5nm')) then nghostcells=max(nghostcells,3) end if if(any(limiter(1:nlevelshi)=='teno5ad')) then nghostcells=max(nghostcells,3) end if if(any(limiter(1:nlevelshi)=='weno5cu6')) then nghostcells=max(nghostcells,3) end if if(any(limiter(1:nlevelshi)=='ppm')) then if(flatsh .or. flatcd) then nghostcells=max(nghostcells,4) else nghostcells=max(nghostcells,3) end if end if if(any(limiter(1:nlevelshi)=='weno7')) then nghostcells=max(nghostcells,4) end if if(any(limiter(1:nlevelshi)=='mpweno7')) then nghostcells=max(nghostcells,4) end if ! If a wider stencil is used, extend the number of ghost cells nghostcells = nghostcells + phys_wider_stencil ! prolongation in AMR for constrained transport MHD needs even number ghosts if(stagger_grid .and. refine_max_level>1 .and. mod(nghostcells,2)/=0) then nghostcells=nghostcells+1 end if select case (coordinate) case (spherical) xprobmin2=xprobmin2*two*dpi;xprobmin3=xprobmin3*two*dpi xprobmax2=xprobmax2*two*dpi;xprobmax3=xprobmax3*two*dpi; case (cylindrical) if (1==phi_) then xprobmin1=xprobmin1*two*dpi;xprobmax1=xprobmax1*two*dpi; end if if (2==phi_) then xprobmin2=xprobmin2*two*dpi;xprobmax2=xprobmax2*two*dpi; end if if (3==phi_) then xprobmin3=xprobmin3*two*dpi;xprobmax3=xprobmax3*two*dpi; end if end select if (log_r0 /= zero) call mpistop("log_r0 offsets the logarithmic radial & &map and is only meaningful for geometry='logSpherical' or & &geometry='logCylindrical'") ! full block size including ghostcells ixGhi1 = block_nx1 + 2*nghostcells ixGhi2 = block_nx2 + 2*nghostcells ixGhi3 = block_nx3 + 2*nghostcells ixGshi1 = ixGhi1 ixGshi2 = ixGhi2 ixGshi3 = ixGhi3 nx_vec = [domain_nx1, domain_nx2, domain_nx3] block_nx_vec = [block_nx1, block_nx2, block_nx3] if (any(nx_vec < 4) .or. any(mod(nx_vec,& 2) == 1)) call mpistop& ('Grid size (domain_nx^D) has to be even and >= 4') if (any(block_nx_vec < 4) .or. any(mod(block_nx_vec,& 2) == 1)) call mpistop& ('Block size (block_nx^D) has to be even and >= 4') if(mod(domain_nx1,block_nx1)/=0) call & mpistop('Grid (domain_nx^D) and block (block_nx^D) must be consistent') if(mod(domain_nx2,block_nx2)/=0) call & mpistop& ('Grid (domain_nx^D) and block (block_nx^D) must be consistent') if(mod(domain_nx3,block_nx3)/=0) call & mpistop& ('Grid (domain_nx^D) and block (block_nx^D) must be consistent') if(refine_max_level>nlevelshi.or.refine_max_level<1)then write(unitterm,*)'Error: refine_max_level',refine_max_level,& '>nlevelshi ',nlevelshi call mpistop("Reset nlevelshi and recompile!") endif dx_vec = [xprobmax1-xprobmin1, xprobmax2-xprobmin2,& xprobmax3-xprobmin3] / nx_vec if (mype==0) then write(c_ndim, '(I1)') 3 write(unitterm, '(A30,' // c_ndim // '(I0," "))') & ' Domain size (cells): ', nx_vec write(unitterm, '(A30,' // c_ndim // '(E9.3," "))') ' Level one dx: ',& dx_vec end if if (any(dx_vec < smalldouble)) call mpistop(& "Incorrect domain size (too small grid spacing)") dx(:, 1) = dx_vec if(sum(w_refine_weight(:))==0) w_refine_weight(1) = 1.d0 if(dabs(sum(w_refine_weight(:))-1.d0)>smalldouble) then write(unitterm,*) "Sum of all elements in w_refine_weight be 1.d0" call mpistop("Reset w_refine_weight so the sum is 1.d0") end if if (mype==0) write(unitterm, '(A30)', advance='no') 'Refine estimation: ' select case (refine_criterion) case (0) if (mype==0) write(unitterm, '(A)') "user defined" case (1) if (mype==0) write(unitterm, '(A)') "relative error" case (2) if (mype==0) write(unitterm, '(A)') "Lohner's original scheme" case (3) if (mype==0) write(unitterm, '(A)') "Lohner's scheme" case default call mpistop("Unknown error estimator, change refine_criterion") end select if (tfixgrid<bigdouble/2.0d0) then if(mype==0)print*,'Warning, at time=',tfixgrid,'the grid will be fixed' end if if (itfixgrid<biginteger/2) then if(mype==0)print*,'Warning, at iteration=',itfixgrid,& 'the grid will be fixed' end if if (ditregrid>1) then if(mype==0)print*,'Note, Grid is reconstructed once every',ditregrid,& 'iterations' end if do islice=1,nslices select case(slicedir(islice)) case(1) ! slicecoord is a physical coordinate; xprobmin1/xprobmax1 hold the ! logical one, which differs from it under LOG_RADIUS. if(slicecoord(islice)<xprobmin1.or.slicecoord(islice)>xprobmax1) & write(uniterr,*)'Warning in read_par_files: ', 'Slice ', islice,& ' coordinate',slicecoord(islice),'out of bounds for dimension ',& slicedir(islice) case(2) if(slicecoord(islice)<xprobmin2.or.slicecoord(islice)>xprobmax2) & write(uniterr,*)'Warning in read_par_files: ', 'Slice ', islice,& ' coordinate',slicecoord(islice),'out of bounds for dimension ',& slicedir(islice) case(3) if(slicecoord(islice)<xprobmin3.or.slicecoord(islice)>xprobmax3) & write(uniterr,*)'Warning in read_par_files: ', 'Slice ', islice,& ' coordinate',slicecoord(islice),'out of bounds for dimension ',& slicedir(islice) end select end do if (mype==0) then write(unitterm, '(A30,A,A)') 'restart_from_file: ', ' ',& trim(restart_from_file) write(unitterm, '(A30,L1)') 'converting: ', convert write(unitterm, '(A)') '' endif deallocate(flux_scheme) end subroutine read_par_files