diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..2f32fc0 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +*.mod +*.o diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..07b018c --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "sac/src/fgdfio"] + path = sac/src/fgdfio + url = git://github.com/SWAT-Sheffield/fgdfio.git diff --git a/sac/par/mhdmodes b/sac/par/binout similarity index 65% rename from sac/par/mhdmodes rename to sac/par/binout index 07ab21b..c9c96b4 100644 --- a/sac/par/mhdmodes +++ b/sac/par/binout @@ -1,18 +1,21 @@ -&testlist / +&testlist + teststr='readfileini' +/ &filelist - filenameini='/archive/mhdmodes_2D.ini' + filenameini='/archive/gdf_testing/mhdmodes_2D.ini' typefileini='binary' - filename= '/archive/mhdmodes_2D.log', - '/archive/mhdmodes_2D.out' + filename= '/archive/gdf_testing/mhdmodes_2D.log', + '/archive/gdf_testing/mhdmodes_2D.out' typefileout='binary' - fullgridout= T + fullgridout= F fullgridini= T / &savelist - dtsave=0.1,0.0001 + dtsave = 1,0.1 + itsave(1,2)=0 / diff --git a/sac/par/binoutmpi b/sac/par/binoutmpi new file mode 100644 index 0000000..12749fa --- /dev/null +++ b/sac/par/binoutmpi @@ -0,0 +1,52 @@ +&testlist + teststr='readfileini' +/ + +&filelist + filenameini='/archive/gdf_testing/mhdmodes_np0202.ini' + + typefileini='binary' + filename= '/archive/gdf_testing/mhdmodes_np0202.log', + '/archive/gdf_testing/mhdmodes_np0202.out' + typefileout='binary' + fullgridout= F + fullgridini= T + / + +&savelist + dtsave = 1,0.1 + itsave(1,2)=0 + + / + + &stoplist + tmax=0.2d0 + itmax = 10000 + / + + &methodlist + + wnames= 'h m1 m2 e b1 b2 eb rhob bg1 bg2' + typefull= 6*'cd4',4*'nul' + typeadvance= 'onestep' + typefilter= 10*'nul' + dimsplit= F + sourcesplit= F + divBfix= F + smallp= 10.d0 + / + + &boundlist + typeB= 10*'fixed' + 10*'fixed' + 10*'fixed' + 10*'fixed' + 10*'fixed' + 10*'fixed' + + / + + ¶mlist + courantpar=0.2 + + / diff --git a/sac/par/gdfini b/sac/par/gdfini new file mode 100644 index 0000000..f106f03 --- /dev/null +++ b/sac/par/gdfini @@ -0,0 +1,52 @@ +&testlist + teststr='readfileini readfileini_gdf' +/ + +&filelist + filenameini='/archive/gdf_testing/mhdmodes_2D_00000000.gdf' + + typefileini='gdf' + filename= '/archive/gdf_testing/mhdmodes_2D.log', + '/archive/gdf_testing/mhdmodes_2D_ii_' + typefileout='gdf' + fullgridout= F + fullgridini= F + / + +&savelist + dtsave=1,0.1 + itsave(1,2)=0 + + / + + &stoplist + tmax=0.2d0 + itmax = 10000 + / + + &methodlist + + wnames= 'h m1 m2 e b1 b2 eb rhob bg1 bg2' + typefull= 6*'cd4',4*'nul' + typeadvance= 'onestep' + typefilter= 10*'nul' + dimsplit= F + sourcesplit= F + divBfix= F + smallp= 10.d0 + / + + &boundlist + typeB= 10*'fixed' + 10*'fixed' + 10*'fixed' + 10*'fixed' + 10*'fixed' + 10*'fixed' + + / + + ¶mlist + courantpar=0.2 + + / diff --git a/sac/par/gdfout b/sac/par/gdfout new file mode 100644 index 0000000..855098a --- /dev/null +++ b/sac/par/gdfout @@ -0,0 +1,52 @@ +&testlist + teststr='readfileini' +/ + +&filelist + filenameini='/archive/gdf_testing/mhdmodes_2D.ini' + + typefileini='binary' + filename= '/archive/gdf_testing/mhdmodes_2D.log', + '/archive/gdf_testing/mhdmodes_2D_' + typefileout='gdf' + fullgridout= F + fullgridini= T + / + +&savelist + dtsave=1,0.1 + itsave(1,2)=0 + + / + + &stoplist + tmax=0.2d0 + itmax = 10000 + / + + &methodlist + + wnames= 'h m1 m2 e b1 b2 eb rhob bg1 bg2' + typefull= 6*'cd4',4*'nul' + typeadvance= 'onestep' + typefilter= 10*'nul' + dimsplit= F + sourcesplit= F + divBfix= F + smallp= 10.d0 + / + + &boundlist + typeB= 10*'fixed' + 10*'fixed' + 10*'fixed' + 10*'fixed' + 10*'fixed' + 10*'fixed' + + / + + ¶mlist + courantpar=0.2 + + / diff --git a/sac/par/gdfoutmpi b/sac/par/gdfoutmpi new file mode 100644 index 0000000..62cc29f --- /dev/null +++ b/sac/par/gdfoutmpi @@ -0,0 +1,52 @@ +&testlist + teststr='' +/ + +&filelist + filenameini='/archive/gdf_testing/mhdmodes_np0202.ini' + + typefileini='binary' + filename= '/archive/gdf_testing/mhdmodes_np0202.log', + '/archive/gdf_testing/mhdmodes_np0202_' + typefileout='gdf' + fullgridout= F + fullgridini= T + / + +&savelist + dtsave = 1,0.1 + itsave(1,2)=0 + + / + + &stoplist + tmax=0.2d0 + itmax = 10000 + / + + &methodlist + + wnames= 'h m1 m2 e b1 b2 eb rhob bg1 bg2' + typefull= 6*'cd4',4*'nul' + typeadvance= 'onestep' + typefilter= 10*'nul' + dimsplit= F + sourcesplit= F + divBfix= F + smallp= 10.d0 + / + + &boundlist + typeB= 10*'fixed' + 10*'fixed' + 10*'fixed' + 10*'fixed' + 10*'fixed' + 10*'fixed' + + / + + ¶mlist + courantpar=0.2 + + / diff --git a/sac/src/Makefile b/sac/src/Makefile index c6f28cf..81fad66 100644 --- a/sac/src/Makefile +++ b/sac/src/Makefile @@ -1,8 +1,7 @@ # This is a VERY striped down makefile for compiling sac with gfortran - -# Check to see if the ENV variables are defined for the compiler +# This will become a SAC only f2py makefile hopefully ifndef F90 - F90=mpif90 + F90=h5pfc endif ifndef F90FLG F90FLG=-std=f2008 @@ -23,11 +22,13 @@ F=.f90 O=.o VACDIR=. +GDFDIR=./fgdfio/lib ################# Definitions for source files ############################# -LIBS = vacdef$F -INCLUDES = vacdef$F vacpar$F vacusrpar$F +LIBS = vacdef$F sacgdf$F +INCLUDES = vacdef$F vacpar$F vacusrpar$F sacgdf$F +GDF_INCLUDES = $(GDFDIR)/gdf_types.F90 $(GDFDIR)/helpers_hdf5.F90 $(GDFDIR)/grid_data_format.F90 $(GDFDIR)/gdf_datasets.F90 VAC_FOR = vac$F vacio$F vacgrid$F vacphys0$F vacphys$F vacusr$F VAC_OBJ = vac$O vacio$O vacgrid$O vacphys0$O vacphys$O vacusr$O @@ -35,8 +36,13 @@ VAC_OBJ = vac$O vacio$O vacgrid$O vacphys0$O vacphys$O vacusr$O # The VACCD, VACMC, VACFCT, VACTVDLF, VACTVD, VACIMPL, VACPOISSON, VACMPI # modules are removed or added in the following two lines by "setvac". # DO NOT TOUCH THESE TWO LINES: -VACFOR = $(VAC_FOR) vaccd$F -VACOBJ = $(VAC_OBJ) vaccd$O +VACFOR = $(VAC_FOR) vaccd$F vacmpi$F +VACOBJ = $(VAC_OBJ) vaccd$O vacmpi$O + +# The VACMPI module is removed or added the following two lines by "setvac". +# DO NOT TOUCH THESE TWO LINES: +VACINIFOR = vacini$F vacio$F vacphys0$F vacusr$F vacmpi$F +VACINIOBJ = vacini$O vacio$O vacphys0$O vacusr$O vacmpi$O ROEOBJ = roetest$O vacphys0$O vacphys$O @@ -47,10 +53,11 @@ PREPROC= $(VACPP_) .SUFFIXES: .SUFFIXES: .t $F $O -$(VACFOR) vacini$F vacdef$F vacpar$F vacusrpar$F : $(PREPROC) +$(VACFOR) vacini$F vacdef$F vacpar$F vacusrpar$F sacgdf$F : $(PREPROC) roetest$F : $(PREPROC) -$(VACOBJ) vacini$O vacall$O vaciniall$O vacsmall$O roetest$O : $(INCLUDES) +$(VACOBJ) vacini$O vacall$O vaciniall$O vacsmall$O roetest$O sacgdf$0 : $(INCLUDES) +$(VACINIOBJ) vacini$O vacall$O vaciniall$O vacsmall$O roetest$O sacgdf$0 : $(INCLUDES) vacusrpar$F: vacusrpar.t $(VACPP) $< $(PREFOR) > $@ @@ -62,6 +69,9 @@ vacpar$F: vacpar.t vacdef$F: vacdef.t $(VACPP) $< $(PREFOR) > $@ +sacgdf$F: sacgdf.t + $(VACPP) $< $(PREFOR) > $@ + $(FOR) -c sacgdf.f90 # General precompilation rule .t$(F): @@ -71,8 +81,6 @@ vacdef$F: vacdef.t $(F)$(O): $(FOR) $(FORFLG) -c $< - - ########### Extra dependencies for some files with "INCLUDE:" statements vacphys$F: vacphys.t vacphys.mhd0.t vacphys.mhdres.t $(VACPP) $< $(PREFOR) > $@ @@ -80,13 +88,20 @@ vacphys$F: vacphys.t vacphys.mhd0.t vacphys.mhdres.t vacgrid$F: vacgrid.t $(VACPP) $< $(PREFOR) > $@ -vac :$(VACOBJ) $(LIBS_) - $(FOR) $(FORFLG) -o $(VACDIR)/vac $(VACOBJ) $(LIBS) +gdf : + $(FOR) $(FORFLG) -c $(GDF_INCLUDES) + +# The "vacini" and "vac" targets use the default Fortran compiler +vacini : gdf $(VACINIOBJ) $(LIBS_) + $(FOR) $(FORFLG) -o $(VACDIR)/vacini $(GDF_INCLUDES) $(VACINIOBJ) $(LIBS) + +vac : gdf $(VACOBJ) $(LIBS_) + $(FOR) $(FORFLG) -o $(VACDIR)/vac $(GDF_INCLUDES) $(VACOBJ) $(LIBS) ###### Removing object files, precompiled Fortran files, and symbolic links clean : - rm -f vac*$F roetest$F Vvac*$F vac*.d *$O *~ + rm -f vac*$F sac*$F roetest$F Vvac*$F vac*.d *$O *~ *.mod cleanall: clean rm -f vacphys.t vacphys0.t vacpar.t vacusrpar.t vacusr.t diff --git a/sac/src/fgdfio b/sac/src/fgdfio new file mode 160000 index 0000000..c66a789 --- /dev/null +++ b/sac/src/fgdfio @@ -0,0 +1 @@ +Subproject commit c66a789eacd9d65beae8d7cf217b35c351ca111f diff --git a/sac/src/sacgdf.t b/sac/src/sacgdf.t new file mode 100644 index 0000000..3c92e9a --- /dev/null +++ b/sac/src/sacgdf.t @@ -0,0 +1,487 @@ +module sacgdf + use hdf5 + +contains + + + subroutine sacgdf_write_file(file_id, gdf_rd, gdf_sp, field_types) + use hdf5, only: HID_T + use gdf, only: gdf_write_file, gdf_root_datasets_T, gdf_parameters_T, gdf_field_type_T + implicit none + + integer(HID_T), intent(inout) :: file_id + type(gdf_root_datasets_T), intent(inout) :: gdf_rd + type(gdf_parameters_T), intent(inout) :: gdf_sp + type(gdf_field_type_T), dimension(:), intent(inout) :: field_types + + call gdf_write_file(file_id, "Sheffield Advanced Code", "GDF Testing version", & + gdf_rd, gdf_sp, field_types) + + call sacgdf_write_eqpar(file_id) + + end subroutine sacgdf_write_file + + + subroutine sacgdf_read_file(file_id, software_name, software_version, & + gdf_rd, gdf_sp, field_types) + use hdf5, only: HID_T + use gdf, only: gdf_read_file, gdf_root_datasets_T, gdf_parameters_T, gdf_field_type_T + implicit none + + integer(HID_T), intent(inout) :: file_id + type(gdf_root_datasets_T), intent(inout) :: gdf_rd + type(gdf_parameters_T), intent(inout) :: gdf_sp + type(gdf_field_type_T), dimension(:), allocatable, intent(inout) :: field_types + character(len=*), intent(out) :: software_name, software_version + + call gdf_read_file(file_id, software_name, software_version, & + gdf_rd, gdf_sp, field_types) + + call sacgdf_read_eqpar(file_id) + + end subroutine sacgdf_read_file + + + subroutine sacgdf_write_eqpar(file_id) + ! Convert simulation parameters to the eqpar array + use common_variables + use hdf5, only: HID_T, h5gopen_f, h5gclose_f + use helpers_hdf5, only: create_attribute + + implicit none + + integer(HID_T), intent(in) :: file_id + + integer(HID_T) :: g_id + integer :: error + real(kind=8), dimension(:), pointer :: r_ptr + integer(kind=4), dimension(:), pointer :: i4_ptr + integer(kind=4), dimension(1), target :: it_arr + + call h5gopen_f(file_id, 'simulation_parameters', g_id, error) + + r_ptr => eqpar(gamma_:gamma_) + call create_attribute(g_id, 'gamma', r_ptr) + + r_ptr => eqpar(eta_:eta_) + call create_attribute(g_id, 'eta', r_ptr) + + r_ptr => eqpar(grav0_:grav0_) + call create_attribute(g_id, 'gravity0', r_ptr) + + r_ptr => eqpar(grav1_:grav1_) + call create_attribute(g_id, 'gravity1', r_ptr) + + ! Read the extra parameters only if we are 2D or 3D + {^IFTWOD + r_ptr => eqpar(grav2_:grav2_) + call create_attribute(g_id, 'gravity2', r_ptr) + + } + {^IFTHREED + r_ptr => eqpar(grav2_:grav2_) + call create_attribute(g_id, 'gravity2', r_ptr + r_ptr => eqpar(grav3_:grav3_) + call create_attribute(g_id, 'gravity3', + } + + r_ptr => eqpar(nu_:nu_) + call create_attribute(g_id, 'nu', r_ptr) + + it_arr(1) = it + i4_ptr => it_arr + call create_attribute(g_id, 'current_iteration', i4_ptr) + + call h5gclose_f(g_id, error) + + end subroutine sacgdf_write_eqpar + + subroutine sacgdf_make_field_types(field_types) + + use common_variables, only: nw + use gdf, only: gdf_field_type_T + implicit none + + type(gdf_field_type_T), dimension(nw), intent(inout) :: field_types + + select case (nw) + case (7) + {^IFONED call sacgdf_make_field_types_1D(field_types) } + case (10) + {^IFTWOD call sacgdf_make_field_types_2D(field_types) } + case (13) + {^IFTHREED call sacgdf_make_field_types_3D(field_types) } + case default + call die("Wrong nw") + end select + + end subroutine sacgdf_make_field_types + + + subroutine sacgdf_make_field_types_1D(field_types) + + use gdf, only: gdf_field_type_T + implicit none + + type(gdf_field_type_T), dimension(7), intent(inout) :: field_types + + ! Write Field Type information + call field_types(1)%init() + field_types(1)%field_to_cgs = 0.001 + field_types(1)%staggering = 0 + field_types(1)%field_units = "kg m^{-3}" + field_types(1)%variable_name = "density_pert" + field_types(1)%field_name = "Pertubation Density" + + call field_types(2)%init() + field_types(2)%field_to_cgs = 0.001 + field_types(2)%staggering = 0 + field_types(2)%field_units = "kg m^{-3}" + field_types(2)%variable_name = "density_bg" + field_types(2)%field_name = "Background Density" + + call field_types(3)%init() + field_types(3)%field_to_cgs = 100 + field_types(3)%staggering = 0 + field_types(3)%field_units = "ms^{-1}" + field_types(3)%variable_name = "velocity_x" + field_types(3)%field_name = "x Component of Velocity" + + call field_types(4)%init() + field_types(4)%field_to_cgs = 10 + field_types(4)%staggering = 0 + field_types(4)%field_units = "Pa" + field_types(4)%variable_name = "internal_energy_pert" + field_types(4)%field_name = "Pertubation Internal Energy" + + call field_types(5)%init() + field_types(5)%field_to_cgs = 10 + field_types(5)%staggering = 0 + field_types(5)%field_units = "Pa" + field_types(5)%variable_name = "internal_energy_bg" + field_types(5)%field_name = "Background Internal Energy" + + call field_types(6)%init() + field_types(6)%field_to_cgs = 10000.0 + field_types(6)%staggering = 0 + field_types(6)%field_units = "T" + field_types(6)%variable_name = "mag_field_x_bg" + field_types(6)%field_name = "x Component of Background Magnetic Field" + + call field_types(7)%init() + field_types(7)%field_to_cgs = 10000.0 + field_types(7)%staggering = 0 + field_types(7)%field_units = "T" + field_types(7)%variable_name = "mag_field_x_pert" + field_types(7)%field_name = "x Component of Pertubation Magnetic Field" + + end subroutine sacgdf_make_field_types_1D + + subroutine sacgdf_make_field_types_2D(field_types) + + use gdf, only: gdf_field_type_T + use common_variables, only: nw + implicit none + + type(gdf_field_type_T), dimension(10), intent(inout) :: field_types + + call sacgdf_make_field_types_1D(field_types(1:7)) + + call field_types(8)%init() + field_types(8)%field_to_cgs = 100 + field_types(8)%staggering = 0 + field_types(8)%field_units = "ms^{-1}" + field_types(8)%variable_name = "velocity_y" + field_types(8)%field_name = "y Component of Velocity" + + call field_types(9)%init() + field_types(9)%field_to_cgs = 10000.0 + field_types(9)%staggering = 0 + field_types(9)%field_units = "T" + field_types(9)%variable_name = "mag_field_y_bg" + field_types(9)%field_name = "y Component of Background Magnetic Field" + + call field_types(10)%init() + field_types(10)%field_to_cgs = 10000.0 + field_types(10)%staggering = 0 + field_types(10)%field_units = "T" + field_types(10)%variable_name = "mag_field_y_pert" + field_types(10)%field_name = "y Component of Pertubation Magnetic Field" + + end subroutine sacgdf_make_field_types_2D + + subroutine sacgdf_make_field_types_3D(field_types) + + use gdf, only: gdf_field_type_T + use common_variables, only: nw + implicit none + + type(gdf_field_type_T), dimension(13), intent(inout) :: field_types + + call sacgdf_make_field_types_2D(field_types(1:10)) + + call field_types(11)%init() + field_types(11)%field_to_cgs = 10000.0 + field_types(11)%staggering = 0 + field_types(11)%field_units = "T" + field_types(11)%variable_name = "mag_field_z_pert" + field_types(11)%field_name = "z Component of Pertubation Magnetic Field" + + call field_types(12)%init() + field_types(12)%field_to_cgs = 10000.0 + field_types(12)%staggering = 0 + field_types(12)%field_units = "T" + field_types(12)%variable_name = "mag_field_z_bg" + field_types(12)%field_name = "z Component of Background Magnetic Field" + + call field_types(13)%init() + field_types(13)%field_to_cgs = 100 + field_types(13)%staggering = 0 + field_types(13)%field_units = "ms^{-1}" + field_types(13)%variable_name = "velocity_z" + field_types(13)%field_name = "z Component of Velocity" + + end subroutine sacgdf_make_field_types_3D + + subroutine sacgdf_read_eqpar(file_id) + ! Convert simulation parameters to the eqpar array + use common_variables + use hdf5, only: HID_T, h5gopen_f, h5gclose_f + use helpers_hdf5, only: read_attribute + + implicit none + + integer(HID_T), intent(in) :: file_id + + integer(HID_T) :: g_id + integer :: error + + real(kind=8), dimension(:), pointer :: r_ptr + integer(kind=4), dimension(:), pointer :: i4_ptr + integer(kind=4), dimension(1), target :: it_arr + + + call h5gopen_f(file_id, 'simulation_parameters', g_id, error) + allocate(r_ptr(1)) + + call read_attribute(g_id, 'gamma', r_ptr) + eqpar(gamma_:gamma_) = r_ptr + + call read_attribute(g_id, 'eta', r_ptr) + eqpar(eta_:eta_) = r_ptr + + call read_attribute(g_id, 'gravity0', r_ptr) + eqpar(grav0_:grav0_) = r_ptr + + call read_attribute(g_id, 'gravity1', r_ptr) + eqpar(grav1_:grav1_) = r_ptr + + ! Read the extra parameters only if we are 2D or 3D + {^IFTWOD + call read_attribute(g_id, 'gravity2', r_ptr) + eqpar(grav2_:grav2_) = r_ptr + } + {^IFTHREED + call read_attribute(g_id, 'gravity2', r_ptr) + eqpar(grav2_:grav2_) = r_ptr + + call read_attribute(g_id, 'gravity3', r_ptr) + eqpar(grav3_:grav3_) = r_ptr + } + + call read_attribute(g_id, 'nu', r_ptr) + eqpar(nu_:nu_) = r_ptr + + allocate(i4_ptr(1)) + call read_attribute(g_id, 'current_iteration', i4_ptr) + it = i4_ptr(1) + + call h5gclose_f(g_id, error) + + end subroutine sacgdf_read_eqpar + + subroutine build_x_array(ix^L, nx, left_edge, right_edge, x) + ! Construct the x array (which is cell centred) from the left_edge and right_edge + ! values, which are also assumed to be cell centred. + implicit none + + integer, intent(in) :: ix^L + integer, dimension(^ND), intent(in) :: nx + real(kind=8), dimension(^ND), intent(in) :: left_edge, right_edge + real(kind=8), dimension(:^D&,:), intent(inout) :: x + + integer :: ix^D + real(kind=8) :: dx^D + +!!$ ! Set ixmin = 1 +!!$ ixmin^D=1; +!!$ ! Set ixmax to ixmin+nx +!!$ ixmax^D=ixmin^D+nx(^D)-1; + + dx^D=(right_edge(^D)-left_edge(^D))/nx(^D); + + { + forall(ix^DD=ixmin^DD:ixmax^DD) + x(ix^DD,^D)= ((ix^D-ixmin^D)*right_edge(^D) + (ixmax^D-ix^D)*left_edge(^D)) / (ixmax^D-ixmin^D) + end forall + \} + + end subroutine build_x_array + + subroutine sacgdf_read_datasets(place, plist_id, w, ix^L) + use hdf5, only: HID_T + use gdf_datasets + use common_variables, only: ixGhi^D, ixGlo^D, nw + use phys_constants, only: mu0 + + implicit none + + integer(HID_T), intent(inout) :: place + integer(HID_T), intent(inout) :: plist_id !< Property list identifier + integer, intent(in) :: ix^L + real(kind=8), intent(inout):: w(ixG^T,nw) + + + call sacgdf_read_datasets_1D(place, plist_id, w, ix^L) + {^IFTWOD call sacgdf_read_datasets_2D(place, plist_id, w, ix^L) } + {^IFTHREED call sacgdf_read_datasets_3D(place, plist_id, w, ix^L) } + + end subroutine sacgdf_read_datasets + + subroutine sacgdf_read_datasets_1D(place, plist_id, w, ix^L) + use hdf5, only: HID_T + use gdf_datasets, only: read_dataset + use common_variables, only: rho_, rhob_, e_, eb_, m1_, b1_, bg1_ + use common_variables, only: ixGhi^D, ixGlo^D, nw, nx + use phys_constants, only: mu0 + + implicit none + + integer(HID_T), intent(inout) :: place + integer(HID_T), intent(inout) :: plist_id !< Property list identifier + integer, intent(in) :: ix^L + real(kind=8), intent(inout) :: w(ixG^T,nw) + + real(kind=8), dimension(:, :, :), allocatable :: wdata3D + + ! Density pert + call read_dataset(place, 'density_pert', wdata3D) + w(ix^S, rho_) = reshape(wdata3D, nx) + deallocate(wdata3D) + + ! Denisty bg + call read_dataset(place, 'density_bg', wdata3D) + w(ix^S, rhob_) = reshape(wdata3D, nx) + deallocate(wdata3D) + + ! Velocity + call read_dataset(place, 'velocity_x', wdata3D) + w(ix^S, m1_) = reshape(wdata3D, nx) * (w(ix^S, rho_) + w(ix^S, rhob_)) + deallocate(wdata3D) + + ! internal energy pert + call read_dataset(place, 'internal_energy_pert', wdata3D) + w(ix^S, e_) = reshape(wdata3D, nx) + deallocate(wdata3D) + + ! internal energy bg + call read_dataset(place, 'internal_energy_bg', wdata3D) + w(ix^S, eb_) = reshape(wdata3D, nx) + deallocate(wdata3D) + + ! Mag field pert + call read_dataset(place, 'mag_field_x_pert', wdata3D) + w(ix^S, b1_) = reshape(wdata3D, nx)/sqrt(mu0) + deallocate(wdata3D) + + ! Mag field bg + call read_dataset(place, 'mag_field_x_bg', wdata3D) + w(ix^S, bg1_) = reshape(wdata3D, nx)/sqrt(mu0) + deallocate(wdata3D) + + end subroutine sacgdf_read_datasets_1D + +{^IFTWOD + subroutine sacgdf_read_datasets_2D(place, plist_id, w, ix^L) + use hdf5, only: HID_T + use gdf_datasets, only: read_dataset + use common_variables, only: rho_, rhob_, m2_, b2_, bg2_ + use common_variables, only: ixGhi^D, ixGlo^D, nw, nx + use phys_constants, only: mu0 + + implicit none + + integer(HID_T), intent(inout) :: place + integer(HID_T), intent(inout) :: plist_id !< Property list identifier + integer, intent(in) :: ix^L + real(kind=8), intent(inout) :: w(ixG^T,nw) + + real(kind=8), dimension(:, :, :), allocatable :: wdata3D + + ! Velocity + call read_dataset(place, 'velocity_y', wdata3D) + w(ix^S, m2_) = reshape(wdata3D, nx) * (w(ix^S, rho_) + w(ix^S, rhob_)) + deallocate(wdata3D) + + ! Mag field pert + call read_dataset(place, 'mag_field_y_pert', wdata3D) + w(ix^S, b2_) = reshape(wdata3D, nx)/sqrt(mu0) + deallocate(wdata3D) + + ! Mag field bg + call read_dataset(place, 'mag_field_y_bg', wdata3D) + w(ix^S, bg2_) = reshape(wdata3D, nx)/sqrt(mu0) + deallocate(wdata3D) + + end subroutine sacgdf_read_datasets_2D +} +{^IFTHREED + subroutine sacgdf_read_datasets_3D(place, plist_id, w, ix^L) + use hdf5, only: HID_T + use gdf_datasets, only: read_dataset + use common_variables, only: rho_, rhob_, m3_, b3_, bg3_ + use common_variables, only: ixGhi^D, ixGlo^D, nw, nx + use phys_constants, only: mu0 + + implicit none + + integer(HID_T), intent(inout) :: place + integer(HID_T), intent(inout) :: plist_id !< Property list identifier + integer, intent(in) :: ix^L + real(kind=8), intent(inout) :: w(ixG^T,nw) + + real(kind=8), dimension(:, :, :), allocatable :: wdata3D + + ! Velocity + call read_dataset(place, 'velocity_y', wdata3D) + w(ix^S, m2_) = reshape(wdata3D, nx) * (w(ix^S, rho_) + w(ix^S, rhob_)) + deallocate(wdata3D) + + ! Mag field pert + call read_dataset(place, 'mag_field_y_pert', wdata3D) + w(ix^S, b2_) = reshape(wdata3D, nx)/sqrt(mu0) + deallocate(wdata3D) + + ! Mag field bg + call read_dataset(place, 'mag_field_y_bg', wdata3D) + w(ix^S, bg2_) = reshape(wdata3D, nx)/sqrt(mu0) + deallocate(wdata3D) + + ! Velocity + call read_dataset(place, 'velocity_z', wdata3D) + w(ix^S, m3_) = reshape(wdata3D, nx) * (w(ix^S, rho_) + w(ix^S, rhob_)) + deallocate(wdata3D) + + ! Mag field pert + call read_dataset(place, 'mag_field_z_pert', wdata3D) + w(ix^S, b3_) = reshape(wdata3D, nx)/sqrt(mu0) + deallocate(wdata3D) + + ! Mag field bg + call read_dataset(place, 'mag_field_z_bg', wdata3D) + w(ix^S, bg3_) = reshape(wdata3D, nx)/sqrt(mu0) + deallocate(wdata3D) + + end subroutine sacgdf_read_datasets_3D +} +end module sacgdf diff --git a/sac/src/vac.t b/sac/src/vac.t index 1acb0ea..fa6c645 100644 --- a/sac/src/vac.t +++ b/sac/src/vac.t @@ -11,21 +11,21 @@ PROGRAM vac ! Pulled upto FORTRAN 2008 by Stuart Mumford 2013 USE constants - USE common_varibles + USE common_variables INTEGER:: ifile,ierrcode,iw - DOUBLE PRECISION:: w(ixG^T,nw),wnrm2,dtold,time0,time1 + REAL(kind=8):: w(ixG^T,nw),wnrm2,dtold,time0,time1 ! functions LOGICAL:: timetofinish,timetosave - DOUBLE PRECISION:: cputime + REAL(kind=8):: cputime !----------------------------------------------------------------------------- {CALL mpiinit ^IFMPI} verbose=.TRUE. .AND.ipe==0^IFMPI IF(verbose)THEN WRITE(*,'(a)')'VAC 4.52 configured to' - WRITE(*,'(a)')' -d=22 -phi=0 -z=0 -g=104,104 -p=mhd -u=stuart' + WRITE(*,'(a)')' -d=22 -phi=0 -z=0 -g=100,100 -p=mhd -u=default' WRITE(*,'(a)')' -on=cd,rk,mpi' WRITE(*,'(a)')' -off=mc,fct,tvdlf,tvd,impl,poisson,ct,gencoord,resist' {^IFMPI WRITE(*,'(a,i3,a)')'Running on ',npe,' processors'} @@ -166,7 +166,7 @@ END PROGRAM vac SUBROUTINE startup USE constants - USE common_varibles + USE common_variables INTEGER:: ifile,iw,ivector,idim,qnvector !----------------------------------------------------------------------------- @@ -235,10 +235,10 @@ SUBROUTINE advance(iws,w) ! Add split sources and fluxes with unsplit sources USE constants - USE common_varibles + USE common_variables INTEGER:: iws(niw_) - DOUBLE PRECISION:: w(ixG^T,nw), w1(ixG^T,nw) + REAL(kind=8):: w(ixG^T,nw), w1(ixG^T,nw) !----------------------------------------------------------------------------- ! Add split sources berforehand if this is required @@ -299,11 +299,11 @@ SUBROUTINE advance_expl(method,ix^L,iws,w1,w) ! w1 can be ised freely. USE constants - USE common_varibles + USE common_variables CHARACTER(^LENTYPE) :: method INTEGER :: ix^L,iws(niw_) - DOUBLE PRECISION :: w(ixG^T,nw),w1(ixG^T,nw) + REAL(kind=8) :: w(ixG^T,nw),w1(ixG^T,nw) LOGICAL :: firstsweep,lastsweep !----------------------------------------------------------------------------- @@ -377,16 +377,16 @@ SUBROUTINE advect(method,ix^L,iws,idim^LIM,w1,w,firstsweep,lastsweep) ! Depending on typeadvance and implpar call advect1 several times USE constants - USE common_varibles + USE common_variables CHARACTER(^LENTYPE):: method INTEGER:: ix^L,iws(niw_),idim^LIM - DOUBLE PRECISION:: w1(ixG^T,nw),w(ixG^T,nw) + REAL(kind=8):: w1(ixG^T,nw),w(ixG^T,nw) ! For most Runge-Kutta type schemes one more full array is needed ! For classical RK4 another array is needed {^ANDIFRK - DOUBLE PRECISION:: w2(ixG^T,nw),w3(ixG^T,nw) + REAL(kind=8):: w2(ixG^T,nw),w3(ixG^T,nw) } !!!MEMORY Needed for typeadvance='adams2' only @@ -502,11 +502,11 @@ SUBROUTINE advect1(method,qdt,ixI^L,iws,idim^LIM,qtC,wCT,qt,w,firstsweep,lastswe ! getboundaries USE constants - USE common_varibles + USE common_variables CHARACTER(^LENTYPE) :: method INTEGER:: ixI^L,ixO^L,iws(niw_),idim^LIM,idim - DOUBLE PRECISION:: qdt,qtC,qt,wCT(ixG^T,nw),w(ixG^T,nw) + REAL(kind=8):: qdt,qtC,qt,wCT(ixG^T,nw),w(ixG^T,nw) LOGICAL, INTENT(INOUT) :: firstsweep,lastsweep !----------------------------------------------------------------------------- @@ -562,10 +562,10 @@ SUBROUTINE addsource2(qdt,ixII^L,ixOO^L,iws,qtC,wCT,qt,w) ! Add source within ixOO for iws: w=w+qdt*S[wCT] USE constants - USE common_varibles + USE common_variables INTEGER:: ixI^L,ixO^L,ixII^L,ixOO^L,iws(niw_) - DOUBLE PRECISION:: qdt,qtC,qt,wCT(ixG^T,nw),w(ixG^T,nw) + REAL(kind=8):: qdt,qtC,qt,wCT(ixG^T,nw),w(ixG^T,nw) !----------------------------------------------------------------------------- oktest=INDEX(teststr,'addsource')>=1 @@ -597,9 +597,9 @@ LOGICAL FUNCTION timetofinish(time0) ! or residual is small enough. Other conditions may be included. USE constants - USE common_varibles + USE common_variables - DOUBLE PRECISION:: time0, cputime + REAL(kind=8):: time0, cputime LOGICAL:: okfinish !----------------------------------------------------------------------------- @@ -622,7 +622,7 @@ LOGICAL FUNCTION timetosave(ifile) ! Other conditions may be included. USE constants - USE common_varibles + USE common_variables INTEGER:: ifile LOGICAL:: oksave @@ -661,9 +661,9 @@ SUBROUTINE getdt_courant(w,ix^L) ! rotations while the value calculated here does not use a rotation. USE constants - USE common_varibles + USE common_variables - DOUBLE PRECISION:: w(ixG^T,nw),cmax(ixG^T),courantmax,dtold + REAL(kind=8):: w(ixG^T,nw),cmax(ixG^T),courantmax,dtold INTEGER:: ix^L,idim LOGICAL:: new_cmax !----------------------------------------------------------------------------- @@ -711,7 +711,7 @@ SUBROUTINE getdt_courant(w,ix^L) END SUBROUTINE getdt_courant !============================================================================= -DOUBLE PRECISION FUNCTION cputime() +REAL(kind=8) FUNCTION cputime() ! Return cputime in seconds as a double precision number. ! For g77 compiler replace F77_ with F77_ everywhere in this function diff --git a/sac/src/vaccd.t b/sac/src/vaccd.t index 18dda7a..252269c 100644 --- a/sac/src/vaccd.t +++ b/sac/src/vaccd.t @@ -13,13 +13,13 @@ SUBROUTINE centdiff4(qdt,ixI^L,ixO^L,iws,idim^LIM,qtC,wCT,qt,w) ! w is the old value at qt on input and the new value at qt+qdt on output. USE constants - USE common_varibles + USE common_variables - DOUBLE PRECISION:: qdt,qtC,qt,wCT(ixG^T,nw),w(ixG^T,nw) + REAL(kind=8):: qdt,qtC,qt,wCT(ixG^T,nw),w(ixG^T,nw) INTEGER:: ixI^L,ixO^L,iws(niw_),idim^LIM LOGICAL :: transport - DOUBLE PRECISION:: v(ixG^T),f(ixG^T), fb(ixG^T) + REAL(kind=8):: v(ixG^T),f(ixG^T), fb(ixG^T) INTEGER:: iiw,iw,ix^L,idim,idir !----------------------------------------------------------------------------- diff --git a/sac/src/vacdef.t b/sac/src/vacdef.t index 962839f..487d67c 100644 --- a/sac/src/vacdef.t +++ b/sac/src/vacdef.t @@ -3,9 +3,10 @@ MODULE phys_constants SAVE INTEGER, PARAMETER :: biginteger=10000000 - DOUBLE PRECISION, PARAMETER :: pi= 3.1415926535897932384626433832795 - DOUBLE PRECISION, PARAMETER :: smalldouble=1.d-99, bigdouble=1.d+99 - DOUBLE PRECISION, PARAMETER :: zero=0d0, one=1d0, two=2d0, half=0.5d0, quarter=0.25d0 + REAL(kind=8), PARAMETER :: pi= 3.1415926535897932384626433832795 + real(kind=8), parameter :: mu0=4.d-7*pi + REAL(kind=8), PARAMETER :: smalldouble=1.d-99, bigdouble=1.d+99 + REAL(kind=8), PARAMETER :: zero=0d0, one=1d0, two=2d0, half=0.5d0, quarter=0.25d0 END MODULE phys_constants @@ -23,7 +24,7 @@ MODULE code_constants INTEGER, PARAMETER :: ixGlo^D=1 ! The next line is edited by SETVAC - INTEGER, PARAMETER :: ixGhi1=104,ixGhi2=104,ixGhimin=104,ixGhimax=104 + INTEGER, PARAMETER :: ixGhi1=100,ixGhi2=100,ixGhimin=100,ixGhimax=100 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INTEGER, PARAMETER :: ndim=^ND, ndir=^NC @@ -89,7 +90,7 @@ MODULE constants END MODULE constants -MODULE common_varibles +module common_variables USE constants SAVE @@ -97,8 +98,8 @@ MODULE common_varibles INTEGER:: ipe, ipe^D, npe, npe^D, nxall^D, nxpe^D, ierrmpi INTEGER:: ixPEmin^D, ixPEmax^D LOGICAL:: mpiupperB(ndim),mpilowerB(ndim) - DOUBLE PRECISION:: sendbuffer(nmpibuffer) - DOUBLE PRECISION:: recvbuffer(nmpibuffer,2) + REAL(kind=8):: sendbuffer(nmpibuffer) + REAL(kind=8):: recvbuffer(nmpibuffer,2) } ! Unit for reading input parameters. @@ -109,7 +110,7 @@ MODULE common_varibles ! General temporary arrays, any subroutine call may change them ! except for subroutines which say the opposite in their header - DOUBLE PRECISION:: tmp(ixG^T),tmp2(ixG^T) + REAL(kind=8):: tmp(ixG^T),tmp2(ixG^T) ! Number of errors during calculation INTEGER:: nerror(nerrcode) @@ -120,36 +121,39 @@ MODULE common_varibles !Grid parameters INTEGER:: ixM^L,ixG^L,nx^D,nx(ndim) INTEGER:: dixB^L + + ! Global x array edges + real(kind=8), dimension(3) :: x_left_edge = (/ 0, 0, 0 /), x_right_edge = (/ 1, 1, 1 /) ! x and dx are local for HPF - DOUBLE PRECISION:: x(IXG^T,ndim),dx(IXG^T,ndim) - DOUBLE PRECISION:: volume,dvolume(IXG^T) - DOUBLE PRECISION:: area(IXGLO1:IXGHI1),areaC(IXGLO1:IXGHI1) - DOUBLE PRECISION:: areadx(IXGLO1:IXGHI1),areaside(IXGLO1:IXGHI1) + REAL(kind=8):: x(IXG^T,ndim),dx(IXG^T,ndim) + REAL(kind=8):: volume,dvolume(IXG^T) + REAL(kind=8):: area(IXGLO1:IXGHI1),areaC(IXGLO1:IXGHI1) + REAL(kind=8):: areadx(IXGLO1:IXGHI1),areaside(IXGLO1:IXGHI1) ! Variables for generalized coordinates and polargrid LOGICAL:: gencoord, polargrid - {^IFGEN DOUBLE PRECISION:: surfaceC(IXG^T,ndim),normalC(IXG^T,ndim,ndim)} - {^NOGEN DOUBLE PRECISION:: surfaceC(2^D&,ndim), normalC(2^D&,ndim,ndim)} + {^IFGEN REAL(kind=8):: surfaceC(IXG^T,ndim),normalC(IXG^T,ndim,ndim)} + {^NOGEN REAL(kind=8):: surfaceC(2^D&,ndim), normalC(2^D&,ndim,ndim)} !Boundary region parameters - DOUBLE PRECISION:: fixB^D(-dixBlo:dixBhi^D%ixGLO^DD:ixGHI^DD,nw) + REAL(kind=8):: fixB^D(-dixBlo:dixBhi^D%ixGLO^DD:ixGHI^DD,nw) INTEGER:: nB,ixB^LIM(ndim,nhiB),idimB(nhiB),ipairB(nhiB) LOGICAL:: upperB(nhiB),fixedB(nw,nhiB),nofluxB(nw,ndim),extraB CHARACTER(^LENTYPE) :: typeB(nw,nhiB),typeBscalar(nhiB) !Equation and method parameters - DOUBLE PRECISION:: eqpar(neqpar+nspecialpar),procpar(nprocpar) + real(kind=8), target :: eqpar(neqpar+nspecialpar),procpar(nprocpar) ! Time step control parameters - DOUBLE PRECISION:: courantpar,dtpar,dtdiffpar,dtcourant(ndim),dtmrpc + REAL(kind=8):: courantpar,dtpar,dtdiffpar,dtcourant(ndim),dtmrpc LOGICAL:: dtcantgrow INTEGER:: slowsteps ! Parameters for the implicit techniques - {^IFPOISSON DOUBLE PRECISION:: wrk(ixG^T,nwrk) } - {^IFIMPL DOUBLE PRECISION:: work(nwork) } + {^IFPOISSON REAL(kind=8):: wrk(ixG^T,nwrk) } + {^IFIMPL REAL(kind=8):: work(nwork) } INTEGER:: nwimpl,nimpl - DOUBLE PRECISION:: implpar,impldiffpar,implerror,implrelax,impldwlimit + REAL(kind=8):: implpar,impldiffpar,implerror,implrelax,impldwlimit INTEGER:: implrestart,implrestart2,impliter,impliternr,implmrpcpar CHARACTER(^LENTYPE) :: typeimplinit,typeimpliter,typeimplmat LOGICAL:: implconserv,implnewton,implcentered,implnewmat @@ -168,23 +172,23 @@ MODULE common_varibles LOGICAL:: divbfix,divbwave,divbconstrain,angmomfix,compactres,smallfix INTEGER:: idimsplit INTEGER:: nproc(nfile+2) - DOUBLE PRECISION:: entropycoef(nw),constraincoef - DOUBLE PRECISION:: smallp,smallpcoeff,smallrho,smallrhocoeff,vacuumrho - DOUBLE PRECISION:: muscleta1,muscleta2,musclomega,acmcoef(nw),acmexpo + REAL(kind=8):: entropycoef(nw),constraincoef + REAL(kind=8):: smallp,smallpcoeff,smallrho,smallrhocoeff,vacuumrho + REAL(kind=8):: muscleta1,muscleta2,musclomega,acmcoef(nw),acmexpo LOGICAL:: acmnolim, fourthorder INTEGER:: acmwidth !Previous time step and residuals - DOUBLE PRECISION:: wold(ixG^T,nw),residual,residmin,residmax + REAL(kind=8):: wold(ixG^T,nw),residual,residmin,residmax ! Flux storage for flux-CT and flux-CD methods !!! for MHD only !!! - {^IFCT DOUBLE PRECISION:: fstore(ixG^T,ndim) } + {^IFCT REAL(kind=8):: fstore(ixG^T,ndim) } !Time parameters INTEGER:: step,istep,nstep,it,itmin,itmax,nexpl,nnewton,niter,nmatvec - DOUBLE PRECISION:: t,tmax,dt,dtmin,cputimemax + REAL(kind=8):: t,tmax,dt,dtmin,cputimemax LOGICAL:: tmaxexact - DOUBLE PRECISION:: tsave(nsavehi,nfile),tsavelast(nfile),dtsave(nfile) + REAL(kind=8):: tsave(nsavehi,nfile),tsavelast(nfile),dtsave(nfile) INTEGER:: itsave(nsavehi,nfile),itsavelast(nfile),ditsave(nfile) INTEGER:: isavet(nfile),isaveit(nfile) @@ -200,6 +204,6 @@ MODULE common_varibles INTEGER:: ixtest1,ixtest2,ixtest3,iwtest,idimtest,ipetest^IFMPI LOGICAL:: oktest !This is a local variable for all subroutines and functions - DOUBLE PRECISION:: maxviscoef + REAL(kind=8):: maxviscoef -END MODULE common_varibles +end module common_variables diff --git a/sac/src/vacgrid.t b/sac/src/vacgrid.t index a29d57e..9746a10 100644 --- a/sac/src/vacgrid.t +++ b/sac/src/vacgrid.t @@ -19,7 +19,7 @@ SUBROUTINE boundsetup ! typeB(iw,iB) - boundary type string USE constants - USE common_varibles + USE common_variables INTEGER:: ix^L,iB,jB,iw,idim,idm,ixG^LIM(ndim),ixM^LIM(ndim) !----------------------------------------------------------------------------- @@ -168,10 +168,10 @@ SUBROUTINE ensurebound(dix,ixI^L,ixO^L,qt,w) ! Adjust ixI and ixO. Call getboundary if needed. USE constants - USE common_varibles + USE common_variables INTEGER:: dix,ixI^L,ixO^L - DOUBLE PRECISION:: qt,w(ixG^T,nw) + REAL(kind=8):: qt,w(ixG^T,nw) !----------------------------------------------------------------------------- oktest=INDEX(teststr,'ensurebound')>0 @@ -197,13 +197,13 @@ END SUBROUTINE ensurebound SUBROUTINE getboundary(qt,iw^LIM,idim^LIM,w) USE constants - USE common_varibles + USE common_variables INTEGER:: iw^LIM,idim^LIM - DOUBLE PRECISION:: qt,w(ixG^T,1:nw) + REAL(kind=8):: qt,w(ixG^T,1:nw) INTEGER:: ix,ix^D,ixe,ixf,ix^L,ixpair^L,idim,iw,iB INTEGER:: iwv,jdim - DOUBLE PRECISION:: coeffnormal,coefftransv + REAL(kind=8):: coeffnormal,coefftransv LOGICAL:: initialized @@ -522,10 +522,10 @@ SUBROUTINE setnoflux(iw,idim,ix^L,fRC,ixR^L,fLC,ixL^L) ! in a boundary region USE constants - USE common_varibles + USE common_variables INTEGER:: iw,idim,ix^L,ixL^L,ixR^L - DOUBLE PRECISION:: fRC(ixG^T), fLC(ixG^T) + REAL(kind=8):: fRC(ixG^T), fLC(ixG^T) INTEGER:: iB,ixe,ixB^L !----------------------------------------------------------------------------- @@ -569,12 +569,12 @@ SUBROUTINE gridsetup1 ! qx - x with an extended index range for calculation of dx USE constants - USE common_varibles + USE common_variables INTEGER:: ix^L,hx^L,jx^L INTEGER:: ix,ixe,ixf,idim,jdim - DOUBLE PRECISION:: qx(IXG^LL^LADD1:,ndim) - DOUBLE PRECISION:: r(IXGLO1-1:IXGHI1+1),rC(IXGLO1-1:IXGHI1+1) + REAL(kind=8):: qx(IXG^LL^LADD1:,ndim) + REAL(kind=8):: r(IXGLO1-1:IXGHI1+1),rC(IXGLO1-1:IXGHI1+1) !----------------------------------------------------------------------------- @@ -671,11 +671,11 @@ END SUBROUTINE gridsetup1 SUBROUTINE gradient4(realgrad,q,ix^L,idim,gradq) USE constants - USE common_varibles + USE common_variables LOGICAL:: realgrad INTEGER:: ix^L,idim - DOUBLE PRECISION:: q(ixG^T),gradq(ixG^T) + REAL(kind=8):: q(ixG^T),gradq(ixG^T) INTEGER:: kx^L,jx^L,hx^L,gx^L INTEGER:: minx1^D,maxx1^D,k !----------------------------------------------------------------------------- @@ -731,10 +731,10 @@ SUBROUTINE laplace4(q,ix^L,laplaceq) !!! We assume uniform Cartesian grid in slab symmetry for now USE constants - USE common_varibles + USE common_variables INTEGER:: ix^L - DOUBLE PRECISION:: q(ixG^T),laplaceq(ixG^T) + REAL(kind=8):: q(ixG^T),laplaceq(ixG^T) INTEGER:: idim,kx^L,jx^L,hx^L,gx^L !----------------------------------------------------------------------------- diff --git a/sac/src/vacini.t b/sac/src/vacini.t index e69de29..48e8fce 100644 --- a/sac/src/vacini.t +++ b/sac/src/vacini.t @@ -0,0 +1,1014 @@ +!############################################################################## +! module vacini + +!INCLUDE:vacnul.process.t +!============================================================================= +program vacini + +use constants + use common_variables + +integer:: iw,ieqpar,ix^L +character*20 :: typeini +double precision:: w(ixG^T,1:nw),wpar(nw) +logical:: lastiw +!----------------------------------------------------------------------------- +{call mpiinit ^IFMPI} + +verbose=.true. .and.ipe==0^IFMPI +if(verbose)then + write(*,'(a)')'VACINI 4.52 configured to' + write(*,'(a)')' -d=33 -phi=0 -z=0 -g=196,54,54 -p= -u=' + write(*,'(a)')' -on=cd,rk,mpi' + write(*,'(a)')' -off=mc,fct,tvdlf,tvd,impl,poisson,ct,gencoord,resist' + {^IFMPI write(*,'(a,i3,a)')'Running on ',npe,' processors'} +endif + +! Some default values +t=zero; it=0; +typefileout='ascii'; typefileini='auto' +snapshotini=0 +fullgridini=.false.; fullgridout=.false. +gencoord= .false. +! There are no ghost cells in VACINI except when "readmesh" is used. +dixB^L=0; ixMmin^D=ixGlo^D; +! Test cell +ixtest^D=ixMmin^D; +! Read parameters from STDIN +unitpar=unitstdin + +{^IFMPI +! MPI reads from a file +unitpar=unitini-1 +open(unitpar,file='vacini.par',status='old') +} + +if(verbose)write(*,*)'Filename for new initial file:' +read(unitpar,'(a)')filename(fileout_) +{^IFMPI +! Extract and check the directional processor numbers and indexes +! and concat the PE number to the output filename +call mpisetnpeDipeD(filename(fileout_)) +} +if(verbose)write(*,*)'Fileheader:' +read(unitpar,'(a)')fileheadout +if(verbose)write(*,*)'Variable names, e.g. "x y rho m1 m2":' +read(unitpar,'(a)')varnames + +call setheaderstrings + +if(verbose)then + write(*,*)'Select action by typing one of the following words: ' + write(*,*)' test' + write(*,*)' typefileini,snapshotini,read,readmesh,readnext,',& + 'typefileout,write,save' + write(*,*)' domain,grid,sheargrid,shiftgrid,polargrid,',& + 'roundgrid,rotategrid' + write(*,*)' transpose,regrid,stretchgrid,stretchgridcent' + write(*,*)' polarvar,spherevar,cartesian,rotatevar,',& + 'setvar,perturbvar,conserve,primitive,multiply,divide' + write(*,*)' uniform,shocktube,wave,wave1,special' + write(*,*)' eqpar,gencoord' +endif + +do + if(verbose)write(*,*)'Action:' + read(unitpar,'(a)')typeini + if(verbose)write(*,*)'> ',typeini + select case(typeini) + case('verbose') + if(verbose)write(*,*)'Verbose:' + read(unitpar,*)verbose + verbose=verbose.and.ipe==0^IFMPI + case('test') + if(verbose)write(*,*)'Teststring:' + read(unitpar,'(a)')teststr + if(verbose)write(*,*)'ixtest, idimtest and iwtest:' + read(unitpar,*) ixtest^D,idimtest,iwtest + {^IFMPI + call mpiix(ixtest^D,ipetest)} + case('typefileini') + if(verbose)write(*,*)'Type of old initial file: ascii/binary/special' + read(unitpar,'(a)')typefileini + case('snapshotini') + if(verbose)write(*,*)'Number of snapshot to be read:' + read(unitpar,*)snapshotini + case('read','readmesh') + if(verbose)write(*,*)'Filename for old initial file:' + read(unitpar,'(a)')filenameini + {call mpisetnpeDipeD(filenameini) ^IFMPI} + if(typeini=='readmesh')then + if(verbose)write(*,*)& + 'Specify boundary width for old initial file:' + read(unitpar,*) dixB^L + fullgridini=.true. + endif + call readfileini(w) + case('readnext') + snapshotini=snapshotini+1 + call readfileini(w) + case('typefileout') + if(verbose)write(*,*)'Type of new initial file: ascii/binary/special' + read(unitpar,'(a)')typefileout + case('write') + call savefile(fileout_,w) + case('save') + call savefile(fileout_,w) + close(unitini+fileout_) + exit + case('grid') + call ini_grid(.true.,ixM^L) + case('domain') + call ini_grid(.false.,ixM^L) + case('shiftgrid') + {^IFMPI call die('shiftgrid not implemented for MPI')} + call shiftgrid(ixM^L,w) + case('sheargrid') + {^IFMPI call die('sheargrid not implemented for MPI')} + call sheargrid(ixM^L,w) + case('polargrid') + {^NOONED call makepolargrid(ixM^L) + gencoord=.true. + if(.false.)}call die('Polar grid is meaningless in 1D') + case('spheregrid') + {^IFTHREED call spheregrid(ixM^L) + gencoord=.true. + if(.false.)}call die('Spherical grid is meaningful in 3D only') + case('roundgrid') + {^NOONED call roundgrid(ixM^L) + gencoord=.true. + if(.false.)}call die('Round grid is meaningless in 1D') + case('rotategrid') + call rotatevar(ixM^L,1,ndim,x) + gencoord=.true. + case('transpose') + {^IFMPI call die('transposexy not implemented for MPI')} + {^IFTWOD call transposexy(ixM^L,w) + if(.false.)}call die('Transpose is implemented in 2D only.') + case('regrid') + {^IFMPI call die('regrid not implemented for MPI')} + call regrid(ixM^L,w) + case('stretchgrid') + {^IFMPI call die('stretchgrid not yet implemented for MPI')} + call stretchgrid(.true.,ixM^L) + case('stretchgridcent') + {^IFMPI call die('stretchgridcent not yet implemented for MPI')} + call stretchgrid(.false.,ixM^L) + case('polarvar') + {^NOONED call polarvar(ixM^L,w) + if(.false.)}call die('Polar variables are meaningless in 1D') + case('spherevar') + {^IFTHREED call spherevar(ixM^L,w) + if(.false.)}call die('Spherical variables are meaningful '// & + 'in 3D only') + case('cartesian') + {^NOONED call cartesian(ixM^L,w) + if(.false.)} call die('Polar variables are meaningless in 1D') + case('rotatevar') + call rotatevar(ixM^L,1,nw,w) + case('setvar') + if(verbose)write(*,*)'Give ix limits:' + read(unitpar,*) ix^L + {^IFMPI + call mpiixlimits(ix^L)} + do + if(verbose)write(*,*) & + 'Variable index, variable value, lastiw (T/F)?' + read(unitpar,*)iw,wpar(iw),lastiw + w(ix^S,iw)=wpar(iw) + if(lastiw)exit + enddo + case('perturbvar') + call perturbvar(ixM^L,w) + case('multiply') + if(verbose)write(*,*)'Give multiplying factors for each variable:' + read(unitpar,*) wpar(1:nw) + do iw=1,nw + w(ixM^S,iw)=w(ixM^S,iw)*wpar(iw) + end do + case('divide') + if(verbose)write(*,*)'Give dividing factors for each variable:' + read(unitpar,*) wpar(1:nw) + do iw=1,nw + w(ixM^S,iw)=w(ixM^S,iw)/wpar(iw) + end do + case('conserv','conserve') + call conserve(ixM^L,w) + case('primitive') + call primitive(ixM^L,w) + case('uniform') + if(verbose)write(*,*)'Give values for each variable:' + read(unitpar,*) wpar(1:nw) + do iw=1,nw + w(ixM^S,iw)=wpar(iw) + end do + case('shocktube') + call ini_shocktube(ixM^L,w) + case('wave') + call ini_wave(ixM^L,w) + case('wave1') + call wave1(ixM^L,w) + case('special') + call specialini(ixM^L,w) + case('eqpar') + if(verbose)write(*,*)'Equation params:',neqpar+nspecialpar + read(unitpar,*)(eqpar(ieqpar),ieqpar=1,neqpar+nspecialpar) + case('gencoord') + if(verbose)write(*,*)'Generalized coordinates (T/F):' + read(unitpar,*)gencoord + case default + call die('Error in VACIni: no such action') + end select +end do + +{^IFMPI +close(unitpar) +call mpifinalize +} + +end + +!============================================================================= +subroutine ini_grid(coord,ix^L) + +! Setup a uniform grid. When coord is .true., the user provides the coordinates +! for the centers of the grid, otherwise the boundaries of the computational +! domaines, thus the centers start at xmin+dx/2, and end at xmax-dx/2 + +use constants + use common_variables + +logical:: coord +integer:: ix^L,ix^D,idim +double precision:: dx^D,xmax(ndim),xmin(ndim) +!----------------------------------------------------------------------------- + +if(verbose)write(*,'(a,3i6)')'Size of mesh. Max: ',ixGhi^D +read(unitpar,*) nx^D +if(verbose)write(*,'(a,3i6)')'Size of mesh: ',nx^D +if(coord)then + if(verbose)write(*,*)'Coordinates of cell centers at the edges' +else + if(verbose)write(*,*)'Boundaries of the computational domain' +endif +if(verbose)write(*,*)'xmin coordinates:' +read(unitpar,*)(xmin(idim),idim=1,ndim) +if(verbose)write(*,*)'xmax coordinates:' +read(unitpar,*)(xmax(idim),idim=1,ndim) + +! Calculate cell sizes and modify coordinates for 'domain' action +if(coord)then + dx^D=(xmax(^D)-xmin(^D))/(nx^D-1); +else + dx^D=(xmax(^D)-xmin(^D))/nx^D; + {^DLOOP + xmax(^D)=xmax(^D)-dx^D/2 + xmin(^D)=xmin(^D)+dx^D/2 + } +endif + +{^IFMPI +! Distribute global grid onto processor cube +nxall^D=nx^D; +call mpigridsetup +} + +ixmax^D=ixmin^D+nx^D-1; +if(ixmax^D>ixGhi^D|.or.) call die('Error in IniGrid: Too big grid') + +{^IFMPI +! Set coordinate limits for this PE +{^DLOOP +xmin(^D) = xmin(^D) + dx^D*(ixPEmin^D-1) +xmax(^D) = xmin(^D) + dx^D*(nx^D-1) +\} +} + +{forall(ix^DD=ixmin^DD:ixmax^DD) x(ix^DD,^D)= & + ((ix^D-ixmin^D)*xmax(^D)+ & + (ixmax^D-ix^D)*xmin(^D)) /(ixmax^D-ixmin^D) \} + +return +end + +!============================================================================= +{^NOONED +subroutine makepolargrid(ix^L) + +use constants + use common_variables + +integer:: ix^L +double precision:: pi2 +!----------------------------------------------------------------------------- + +if(verbose)write(*,*)& + 'First coordinate is interpreted as radius, second as angle/2pi.' + +pi2=8*atan(one) + +tmp(ix^S)=x(ix^S,1) +x(ix^S,1)=tmp(ix^S)*cos(x(ix^S,2)*pi2) +x(ix^S,2)=tmp(ix^S)*sin(x(ix^S,2)*pi2) + +return +end +} +!============================================================================= +{^IFTHREED +subroutine spheregrid(ix^L) + +use constants + use common_variables + +integer:: ix^L +double precision:: pi2 +!----------------------------------------------------------------------------- + +if(verbose)write(*,*)'R, PHI/2Pi [0,1], THETA/2Pi [-0.25,0.25] --> X,Y,Z' + +pi2=8*atan(one) + +tmp(ix^S)=x(ix^S,1) +x(ix^S,1)=tmp(ix^S)*cos(x(ix^S,3)*pi2)*cos(x(ix^S,2)*pi2) +x(ix^S,2)=tmp(ix^S)*cos(x(ix^S,3)*pi2)*sin(x(ix^S,2)*pi2) +x(ix^S,3)=tmp(ix^S)*sin(x(ix^S,3)*pi2) + +return +end +} +!============================================================================= +{^NOONED +subroutine roundgrid(ix^L) + +! Calculate the shrink factor to shrink a rectangle to an ellipse. For a sguare +! 1 in directions parallel to x and y +! 1-r+r*sqrt(0.5) in diagonal directions, where r is the radial distance +! normalized to 1. + +use constants + use common_variables + +integer:: ix^L +double precision:: dist1(ixG^T),dist2(ixG^T),weight(ixG^T) +double precision:: xcent^D,rounded,squared +!----------------------------------------------------------------------------- + +if(verbose)write(*,*)'Coordinates of center point:' +read(unitpar,*) xcent^D +if(verbose)then + write(*,*)'Center of rounded grid:',xcent^D + write(*,*)'If the rectangle is mapped to the (-1,-1,1,1) square' + write(*,*)'give the distances to be rounded, and to be squared:' +endif +read(unitpar,*)rounded,squared + +! Normalized distances in the L1 and L2 norms + +dist1(ix^S)=max(^D&abs((x(ix^S,^D)-xcent^D)/(x(ixmax^DD,^D)-xcent^D))) +dist2(ix^S)=sqrt(^D&((x(ix^S,^D)-xcent^D)/(x(ixmax^DD,^D)-xcent^D))**2+) + +! The weight increases from 0 to 1 for distance between 0 and rounded, and +! it drops back to 0 in the distance range rounded and squared. +where(dist1(ix^S)nw.or.iw2<1.or.iw2>nw)call die( & + 'Error in RotateVar: Incorrect iw1 and iw2') + +w1(ix^S)=w(ix^S,iw1) +w2(ix^S)=w(ix^S,iw2) +w(ix^S,iw1)=cos(angle)*w1(ix^S)-sin(angle)*w2(ix^S) +w(ix^S,iw2)=sin(angle)*w1(ix^S)+cos(angle)*w2(ix^S) + +return +end + +{^IFTWOD +!============================================================================= +subroutine transposexy(ix^L,w) + +! Transpose the first two coordinates of the grid (x), the variables (w), +! and exchange the vector components as required + +use constants + use common_variables + +integer:: ix^L,ixold^L,ix^D,iw,ivect,idim,qnvector +double precision:: w(ixG^T,1:nw) +!----------------------------------------------------------------------------- +ixold^L=ix^L; +ix^LIM1=ixold^LIM2; +ix^LIM2=ixold^LIM1; +! Transpose x +do idim=1,ndim + !!!For sake of f90tof77 x(ix^S,idim)=transpose(x(ixold^S,idim)) is replaced: + tmp(ix^S)=x(ix^S,idim) + {do ix^D=ixmin^D,ixmax^D\} + x(ix^D,idim)=tmp(ix^DB) + {enddo^D&\} +enddo + +! Swap the X and Y coordinates +tmp(ix^S)=x(ix^S,1) +x(ix^S,1)=x(ix^S,2) +x(ix^S,2)=tmp(ix^S) + +! Transpose w +do iw=1,nw + !!!For sake of f90tof77 w(ix^S,iw)=transpose(w(ixold^S,iw)) is replaced by + tmp(ix^S)=w(ix^S,iw) + {do ix^D=ixmin^D,ixmax^D\} + w(ix^D,iw)=tmp(ix^DB) + {enddo^D&\} +enddo + +! Swap the vector variables +! qnvector is only used to avoid compiler warning when nvector=0 +qnvector=nvector +do ivect=1,qnvector + if(verbose)write(*,"(a,i1,a)")& + 'Index of first component of vector variable #',ivect,':' + read(unitpar,*)iw + if(iw>=nw.or.iw<1)call die('Error in TransposeXY: Incorrect iw.') + tmp(ix^S)=w(ix^S,iw) + w(ix^S,iw)=w(ix^S,iw+1) + w(ix^S,iw+1)=tmp(ix^S) +end do + +return +end +} +!============================================================================= +subroutine regrid(ix^L,w) + +! Change the number of grid points and extrapolate and interpolate the +! original cell positions x and averaged values w. + +use constants + use common_variables + +integer:: ix^L,nix^D,ixnew^L,iw,idim +double precision:: w(ixG^T,1:nw),q(ixG^T) +!----------------------------------------------------------------------------- + +if(verbose)write(*,*)'Define number of gridpoints in each direction:' +read(unitpar,*) nix^D + +ixnewmin^D=ixmin^D; +ixnewmax^D=ixmin^D+nix^D-1; + +do idim=1,ndim + q(ix^S)=x(ix^S,idim) + call regrid1(ix^L,ixnew^L,q) + x(ixnew^S,idim)=q(ixnew^S) +enddo +do iw=1,nw + q(ix^S)=w(ix^S,iw) + call regrid1(ix^L,ixnew^L,q) + w(ixnew^S,iw)=q(ixnew^S) +enddo + +ixmax^D=ixnewmax^D; + +return +end + +!============================================================================= +subroutine regrid1(ixI^L,ixO^L,q) + +! Interpolate q from grid determined by ixI to ixO. Use the distances measured +! between the Cartesian gridpoints, i.e. distances in generalized coordinates. +! The DOMAIN ixImin-0.5..ixImax+0.5 is rediscretized by ixOmax-ixOmin+1 points. + +use constants + use common_variables + +integer:: ixI^L,ixO^L,ixI^D,ixO^D,dixI^D +double precision:: q(ixG^T),qnew(ixG^T),dxO^D,xO^D,coeff^D(0:1) +!----------------------------------------------------------------------------- + +! Grid spacing of the output grid stretched onto the integer input grid +dxO^D=(ixImax^D-ixImin^D+one)/(ixOmax^D-ixOmin^D+one); + +qnew(ixO^S)=zero +{do ixO^D=ixOmin^D,ixOmax^D\} + + ! Location of the output grid point + xO^D=ixImin^D-half+(ixO^D-ixOmin^D+half)*dxO^D; + + ! Index of the input grid point to the left of xO within ixImin..ixImax-1 + ixI^D=min(ixImax^D-1,max(ixImin^D,int(xO^D))); + + ! Calculate bilinear interpolation/extrapolation coefficients + coeff^D(1)=xO^D-ixI^D; + coeff^D(0)=1-coeff^D(1); + + ! Interpolate q into qnew + {do dixI^D=0,1\} + qnew(ixO^D)=qnew(ixO^D)+(coeff^D(dixI^D)*)*q(ixI^D+dixI^D) + {enddo^D&\} + +{enddo^D&\} + +q(ixO^S)=qnew(ixO^S); + +return +end + +!============================================================================= +subroutine stretchgrid(qdomain,ix^L) + +! Stretch the grid logarithmically in direction idim segment by segment +! The original computational domain size is preserved if qdomain is true, +! and the first and last grid center locations are preserved if it is false. + +use constants + use common_variables + +logical:: qdomain +integer:: ix^L,ix,ixL,ixR,idim,iseg,nseg +integer,parameter:: qixhi=10000 +double precision:: qxL,qxR,qdxL,qdxR,qdxsum,qdx(qixhi) +!----------------------------------------------------------------------------- + +if(verbose)write(*,*)'Direction of stretch and number of segments' +read(unitpar,*)idim,nseg + +if(idim>ndim.or.idim<1)call die('Error in StretchGrid: Invalid direction') +if(nseg<1)call die('Error in StretchGrid: Invalid number of segments') + +if(qdomain)then + qxL=1.5*x(ixmin^D,idim)-0.5*x(ixmin^D+1,idim) + qxR=1.5*x(ixmax^D,idim)-0.5*x(ixmax^D-1,idim) + if(verbose)write(*,*)'Old domain from',qxL,' to',qxR +else + qxL=x(ixmin^D,idim) + qxR=x(ixmax^D,idim) + if(verbose)write(*,*)'Old centers from',qxL,' to',qxR +endif + +if(qxL>=qxR)& + call die('Error in StretchGrid: qxRqixhi)call die( & + 'Error in StretchGrid: Too big grid, change qixhi') + ixL=ixmin^D+1 + qdxL=one + qdx(ixL)=qdxL + if(qdomain)then + qdxsum=1.5D0 + else + qdxsum=1.0D0 + endif + do iseg=1,nseg + if(verbose)write(*,*)'xlast-xprev for segment:',iseg + read(unitpar,*)qdxR + if(iseg=ixmax^D) & + call die('Error in StretchGrid: bad segment position') + else + if(verbose)write(*,*)'Last segment ends at right boundary' + ixR=ixmax^D + end if + do ix=ixL+1,ixR + qdx(ix)=qdxL*(qdxR/qdxL)**(dble(ix-ixL)/dble(ixR-ixL)) + qdxsum=qdxsum+qdx(ix) + enddo + qdxL=qdxR + ixL=ixR + end do + if(qdomain)qdxsum=qdxsum+half*qdx(ixR) + qdx(ixmin^D+1:ixmax^D)=qdx(ixmin^D+1:ixmax^D)*(qxR-qxL)/qdxsum + if(qdomain)x(ixmin^D^D%ix^S,idim)=qxL+half*qdx(ixmin^D+1) + do ix=ixmin^D+1,ixmax^D + x(ix^D%ix^S,idim)=x(ix-1^D%ix^S,idim)+qdx(ix) + enddo + \} +end select + +if(qdomain)then + qxL=1.5*x(ixmin^D,idim)-0.5*x(ixmin^D+1,idim) + qxR=1.5*x(ixmax^D,idim)-0.5*x(ixmax^D-1,idim) + if(verbose)write(*,*)'New domain from',qxL,' to',qxR +else + qxL=x(ixmin^D,idim) + qxR=x(ixmax^D,idim) + if(verbose)write(*,*)'New centers from',qxL,' to',qxR +endif + +return +end + +!============================================================================= +subroutine perturbvar(ix^L,w) + +! Perturb a variable within limits ixP by adding the product of sine waves +! in each direction. The phases of the waves are relative to ixPmin. + +use constants + use common_variables + +integer:: ix^L,ixP^L,idim,iw +double precision:: w(ixG^T,nw),dw,wavenum(ndim),phase(ndim) +!----------------------------------------------------------------------------- + +if(verbose)write(*,*)'Variable index, amplitude and',& + ' ixP^DL ' +read(unitpar,*)iw,dw,ixP^DL +{^IFMPI +call mpiixlimits(ixP^L)} +if(verbose)write(*,*)'wavenumber and phase for each idim:' +read(unitpar,*)(wavenum(idim),phase(idim),idim=1,ndim) + +w(ixP^S,iw)=w(ixP^S,iw)+dw* & + {sin(phase(^D)+(x(ixP^S,^D)-x(ixPmin^DD,^D))*wavenum(^D))*} + +return +end + +!============================================================================= +subroutine ini_shocktube(ix^L,w) + +! The shocktube is divided into nseg segments in the chosen idim direction. +! Linear interpolation in segments with given left and right states. + +use constants + use common_variables + +integer:: ix^L +double precision:: w(ixG^T,nw) +integer:: ix,ixL,ixR,idim,iw,iseg,nseg,ieqpar +double precision:: wL(nw),wR(nw) +{^IFMPI logical:: inside} +!----------------------------------------------------------------------------- + +if(verbose)write(*,*)'Normal of slab symmetry' +read(unitpar,*)idim +if(verbose)write(*,*)'Number of segments' +read(unitpar,*)nseg +if(verbose)write(*,*)'All variables at minimal position:' +read(unitpar,*)wL(1:nw) +select case(idim) + {case(^D) + ixL=ixmin^D + {ixL=ixL-ixPEmin^D+1 ^IFMPI} + do iseg=1,nseg + if(verbose)write(*,*)'All variables at end of segment:',iseg + read(unitpar,*)wR(1:nw) + if(iseg=ixmin^D.and.ix<=ixmax^D) & + w(ix^D%ix^S,iw)=((ixR-ix)*wL(iw)+(ix-ixL)*wR(iw))/(ixR-ixL) + + !w(ix^D%ix^S,iw)=((x(ixR^D%ix^S,^D)-x(ix^D%ix^S,^D))*wL(iw)+ & + ! (x(ix^D%ix^S,^D)-x(ixL^D%ix^S,^D))*wR(iw))& + ! /(x(ixR^D%ix^S,^D)-x(ixL^D%ix^S,^D)) + end do + end do + ixL=ixR + wL(1:nw)=wR(1:nw) + end do \} + case default + call die('Error in Ini_ShockTube: Unknown dimension') +end select + +if(verbose)write(*,*)'Eqpar:' +read(unitpar,*)(eqpar(ieqpar),ieqpar=1,neqpar+nspecialpar) + +return +end + +!============================================================================= +subroutine ini_wave(ix^L,w) + +! Sum of sine waves in each direction and each variable + +use constants + use common_variables + +integer:: ix^L,ix^D,idim,iw,ieqpar +double precision:: w(ixG^T,nw),wmean(nw),wavenum(ndim),ampl(ndim),phase(ndim) +logical:: nextiw +!----------------------------------------------------------------------------- + +if(verbose)write(*,*)'Mean values:' +read(unitpar,*)wmean(1:nw) +do iw=1,nw + w(ix^S,iw)=wmean(iw) +end do +do iw=1,nw + if(verbose)write(*,*)'Variable iw:',iw + do + if(verbose)write(*,*)'Amplitude,wavenum,phase for each dim,nextiw=T/F:' + read(unitpar,*)(ampl(idim),wavenum(idim),phase(idim),idim=1,ndim),nextiw + do idim=1,ndim + w(ix^S,iw)=w(ix^S,iw)+ampl(idim)*sin(phase(idim)+& + x(ix^S,idim)*wavenum(idim)) + enddo + if(nextiw)exit + enddo +enddo + +if(verbose)write(*,*)'Eqpar:' +read(unitpar,*)(eqpar(ieqpar),ieqpar=1,neqpar+nspecialpar) + +return +end + +!============================================================================= +subroutine wave1(ix^L,w) + +! Sine waves with arbitrary wave vectors and shifts using rationalized angle +! units (1.0 = full circle) + +use constants + use common_variables + +integer:: ix^L,idim,iw +double precision:: w(ixG^T,nw),wmean,ampl,wavenum(ndim),phase,pi2 +logical:: lastiw +!----------------------------------------------------------------------------- + +pi2=8*atan(one) +if(verbose)& + write(*,*)'w(iw)=wmean+dw*sin(2*pi*[x*kx+y*ky+phase]) (Note the 2*pi!)' +do + if(verbose)write(*,*)'Give iw,wmean,dw,k(idim),phase,lastiw (quit with T):' + read(unitpar,*)iw,wmean,ampl,(wavenum(idim),idim=1,ndim),phase,lastiw + w(ix^S,iw)=wmean+ampl*sin(pi2*({wavenum(^D)*x(ix^S,^D)+}+phase)) + if(lastiw)exit +enddo + +return +end +!============================================================================= +! Some interface routines for subroutines often used in the VACUSR module +! to keep the compiler happy +!============================================================================= +subroutine gradient(realgrad,q,ix^L,idir,gradq) +logical:: realgrad +integer:: ix^L,idir +double precision:: q(*),gradq(*) +call die('Error: VACINI cannot call gradient !') +end + +!============================================================================= +subroutine laplace4(q,ix^L,laplaceq) +integer:: ix^L +double precision:: q(*),laplaceq(*) +call die('Error: VACINI cannot call laplace4 !') +end + +!============================================================================= +subroutine ensurebound(dix,ixI^L,ixO^L,qt,w) +integer:: dix,ixI^L,ixO^L +double precision:: qt,w(*) +call die('Error: VACINI cannot call ensurebound !') +end +!============================================================================= +! end module vacini +!############################################################################## + + diff --git a/sac/src/vacio.t b/sac/src/vacio.t index db17bae..c47e6e9 100644 --- a/sac/src/vacio.t +++ b/sac/src/vacio.t @@ -1,23 +1,25 @@ + !############################################################################## ! module vacio !============================================================================= -SUBROUTINE readparameters(w) +subroutine readparameters(w) ! This subroutine sets or reads all default parameters from par/DEFAULT, ! then it reads the par/PROBLEM parameter file through standard input, ! and the initial data from data/PROBLEM.ini as soon as the filename is read ! from the parameter file. - USE constants - USE common_varibles + use constants + use common_variables + + real(kind=8):: w(ixG^T,nw) - DOUBLE PRECISION:: w(ixG^T,nw) + character(^LENTYPE):: typepred(nw),typefull(nw),typeimpl(nw),typefilter(nw), mpifiletype - CHARACTER(^LENTYPE):: typepred(nw),typefull(nw),typeimpl(nw),typefilter(nw) - DOUBLE PRECISION:: muscleta - INTEGER:: i,j,k,iw,idim,iB,ifile,isave - LOGICAL:: implmrpc,globalixtest^IFMPI + real(kind=8):: muscleta + integer:: i,j,k,iw,idim,iB,ifile,isave + logical:: implmrpc,globalixtest^IFMPI ! The use of NAMELIST is not recommended by Fortran 90. It could be replaced ! by some string manipulations, but that is difficult in Fortran 77/Adaptor. @@ -25,15 +27,15 @@ SUBROUTINE readparameters(w) ! SM: NAMELIST is a defined part of the F90+ specification. - NAMELIST /testlist/ teststr,ixtest1,ixtest2,ixtest3,& + namelist /testlist/ teststr,ixtest1,ixtest2,ixtest3,& iwtest,idimtest,ipetest^IFMPI - NAMELIST /filelist/ filenameini,filename,varnames{,^IFMPI npe^D}, & + namelist /filelist/ filenameini,filename,varnames{,^IFMPI npe^D}, & typefileini,typefileout,typefilelog,& snapshotini,snapshotout,fullgridini,fullgridout,dixB^L - NAMELIST /savelist/ tsave,itsave,dtsave,ditsave - NAMELIST /stoplist/ itmax,tmax,tmaxexact,dtmin,residmin,residmax,t,it,& + namelist /savelist/ tsave,itsave,dtsave,ditsave + namelist /stoplist/ itmax,tmax,tmaxexact,dtmin,residmin,residmax,t,it,& cputimemax - NAMELIST /methodlist/ wnames,fileheadout,eqpar,& + namelist /methodlist/ wnames,fileheadout,eqpar,& typeadvance,typefull,typepred,typeimpl,typefilter,& typelimiter,typeentropy,entropycoef,artcomp,& typelimited,typefct,typetvd,typeaxial,& @@ -46,8 +48,8 @@ SUBROUTINE readparameters(w) divbfix,divbwave,divbconstrain,angmomfix,compactres,& smallfix,smallp,smallpcoeff,smallrho,smallrhocoeff,& vacuumrho,nproc,procpar - NAMELIST /boundlist/ nB,typeB,ixB^LIM,idimB,upperB,extraB,typeBscalar,ipairB - NAMELIST /paramlist/ courantpar,dtpar,dtdiffpar,dtcantgrow,slowsteps,& + namelist /boundlist/ nB,typeB,ixB^LIM,idimB,upperB,extraB,typeBscalar,ipairB + namelist /paramlist/ courantpar,dtpar,dtdiffpar,dtcantgrow,slowsteps,& implmrpcpar,& implpar,impldiffpar,implerror,implrelax,impldwlimit,& implrestart,implrestart2,impliter,impliternr,& @@ -63,191 +65,194 @@ SUBROUTINE readparameters(w) ! Set default values for arrays (except the ones read from the inifile) - DO ifile=1,nfile - DO isave=1,nsavehi + do ifile=1,nfile + do isave=1,nsavehi tsave(isave,ifile)=bigdouble ! t of saves into the output files itsave(isave,ifile)=biginteger ! it of saves into the output files - END DO + end do dtsave(ifile)=bigdouble ! time between saves ditsave(ifile)=biginteger ! timesteps between saves isavet(ifile)=1 ! index for saves by t isaveit(ifile)=1 ! index for saves by it - END DO + end do - DO iw=1,nw + do iw=1,nw typepred(iw)='default' ! Predictor scheme (will be adjusted later) typefull(iw)='tvdlf' ! Full step scheme typeimpl(iw)='nul' ! Implicit step scheme typefilter(iw)='nul' ! Filter scheme typelimiter(iw)='minmod' ! Limiter type for flow variables/characteristics typeentropy(iw)='nul' ! Entropy fix type - artcomp(iw)=.FALSE. ! No artificial compression for Harten type TVD - END DO + artcomp(iw)=.false. ! No artificial compression for Harten type TVD + end do - DO iw=1,nw + do iw=1,nw acmcoef(iw)=-one ! Coefficients (0,1) for the dissipative fluxes - ENDDO ! negative value means no coefficient is used + enddo ! negative value means no coefficient is used - DO i=1,nfile+2 ! Elements define processing for the fullstep, + do i=1,nfile+2 ! Elements define processing for the fullstep, nproc(i)=0 ! halfstep and the nfile output files. If the value - END DO ! is 0, no processing. For nproc(1) and nproc(2) + end do ! is 0, no processing. For nproc(1) and nproc(2) ! the value defines the proc. frequency. Negative ! value results in a call at every sweep. Positive ! value N results in a call at every N-th step before ! the first sweep. For nproc(ifile+2) the nonzero ! values cause processing for that file. - DO i=1,nprocpar + do i=1,nprocpar procpar(i)=-one ! Parameters for processing - END DO + end do smallp=-one ! Default small pressure, redefined in getpthermal smallrho=-one ! Default small density, redefined in keeppositive vacuumrho=-one ! Density representing vacuum nB=2*ndim ! If nB is not specified by the user, gridsetup - DO iB=1,nhiB ! will create 2*ndim boundary regions by default. - DO iw=1,nw + do iB=1,nhiB ! will create 2*ndim boundary regions by default. + do iw=1,nw typeB(iw,iB) ='cont' ! Default boundary type - fixedB(iw,iB)=.FALSE. ! Fixed boundaries are not extrapolated into yet - END DO + fixedB(iw,iB)=.false. ! Fixed boundaries are not extrapolated into yet + end do ipairB(iB)=0 ! periodic pair is unknown, but can be set or guessed - END DO - DO iw=1,nw - DO idim=1,ndim - nofluxB(iw,idim)=.FALSE. ! No zero flux condition for variables - ENDDO - END DO + end do + do iw=1,nw + do idim=1,ndim + nofluxB(iw,idim)=.false. ! No zero flux condition for variables + enddo + end do dixB^L=2; ! Default width of boundary regions ixBmax(1,1)=0 ! An impossible value if user specifies boundaries. ! Read scalar parameters from the par/DEFAULT file unitpar=unitini-1 - OPEN(unitpar,file='par/DEFAULT',status='old') + open(unitpar,file='par/DEFAULT',status='old') - READ(unitpar,testlist) - READ(unitpar,filelist) - READ(unitpar,savelist) - READ(unitpar,stoplist) - READ(unitpar,methodlist) - READ(unitpar,boundlist) - READ(unitpar,paramlist) + read(unitpar,testlist) + read(unitpar,filelist) + read(unitpar,savelist) + read(unitpar,stoplist) + read(unitpar,methodlist) + read(unitpar,boundlist) + read(unitpar,paramlist) - CLOSE(unitpar) + close(unitpar) ! end defaults ! Initialize Kronecker delta, and Levi-Civita tensor - DO i=1,3 - DO j=1,3 - IF(i==j)THEN + do i=1,3 + do j=1,3 + if(i==j)then kr(i,j)=1 - ELSE + else kr(i,j)=0 - ENDIF - DO k=1,3 - IF(i==j.OR.j==k.OR.k==i)THEN + 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 + else if(i+1==j.or.i-2==j)then lvc(i,j,k)=1 - ELSE + else lvc(i,j,k)=-1 - ENDIF - ENDDO - ENDDO - ENDDO + endif + enddo + enddo + enddo ! Initialize error conunters and equation parameters - DO i=1,nerrcode + do i=1,nerrcode nerror(i)=0 - END DO - DO i=1,neqpar+nspecialpar + end do + do i=1,neqpar+nspecialpar eqpar(i)=zero - END DO + end do ! read from STDIN unitpar=unitstdin {^IFMPI ! MPI reads from a file unitpar=unitini-1 - OPEN(unitpar,file='vac.par',status='old') + open(unitpar,file='vac.par',status='old') } ! Start reading parameters from standard input, i.e. "< par/PROBLEM" {ipetest = -1 ^IFMPI} - READ(unitpar,testlist) + read(unitpar,testlist) {^IFMPI ! ixtest^D is given for the full grid if ipetest was not set explicitly globalixtest = ipetest < 0 - ipetest = MAX(ipetest,0) + ipetest = max(ipetest,0) ! Erase test string for other processors unless ipetest >= npe (test all PEs) - IF(ipetest=1 - IF(oktest) WRITE(unitterm,*)'ReadParameters' - IF(oktest) WRITE(unitterm,testlist) + oktest=index(teststr,'readparameters')>=1 + if(oktest) write(unitterm,*)'ReadParameters' + if(oktest) write(unitterm,testlist) varnames='default' - READ(unitpar,filelist) + read(unitpar,filelist) {^IFMPI ! Extract and check the directional processor numbers and indexes ! and concat the PE number to the input and output filenames - CALL mpisetnpeDipeD(filenameini) - CALL mpisetnpeDipeD(filename(fileout_)) + mpifiletype = 'inifile' + call mpisetnpeDipeD(filenameini, mpifiletype) + mpifiletype = 'outfile' + call mpisetnpeDipeD(filename(fileout_), mpifiletype) } - IF(oktest) THEN - {^IFMPI WRITE(unitterm,*)'npe^D=',npe^D} - WRITE(unitterm,*)filenameini - DO ifile=1,nfile - WRITE(unitterm,*)filename(ifile) - ENDDO - WRITE(unitterm,*)'Type of ini/out and log files:',& + if(oktest) then + {^IFMPI write(unitterm,*)'npe^D=',npe^D} + write(unitterm,*)filenameini + do ifile=1,nfile + write(unitterm,*)filename(ifile) + enddo + write(unitterm,*)'Type of ini/out and log files:',& typefileini,typefileout,typefilelog - IF(varnames/='default')WRITE(unitterm,*)'Varnames:',varnames - IF(snapshotini>0)WRITE(unitterm,*)'Snapshotini:',snapshotini - IF(snapshotout>0)WRITE(unitterm,*)'Snapshotout:',snapshotout - WRITE(unitterm,*)'Fullgridini,out:',fullgridini,fullgridout - ENDIF + if(varnames/='default')write(unitterm,*)'Varnames:',varnames + if(snapshotini>0)write(unitterm,*)'Snapshotini:',snapshotini + if(snapshotout>0)write(unitterm,*)'Snapshotout:',snapshotout + write(unitterm,*)'Fullgridini,out:',fullgridini,fullgridout + endif - CALL readfileini(w) + call readfileini(w) + print*, "back in readparams" ! Default for output header line fileheadout=fileheadini {^IFMPI ! Reset global test cell indexes to local ones - IF(globalixtest)CALL mpiix(ixtest^D,ipetest) + if(globalixtest)call mpiix(ixtest^D,ipetest) } - READ(unitpar,savelist) - DO ifile=1,nfile - IF(dtsave(ifile)=1 + oktest=index(teststr,'readfileini')>=1 - IF(oktest) WRITE(unitterm,*)'ReadFileIni' + if(oktest) write(unitterm,*)'ReadFileIni' - INQUIRE(file=filenameini,exist=fileexist) - IF(.NOT.fileexist) CALL die('Stop: file does not exist, filenameini='//& + inquire(file=filenameini,exist=fileexist) + if(.not.fileexist) call die('Stop: file does not exist, filenameini='//& filenameini) - OPEN(unitini,file=filenameini,status='old') + open(unitini,file=filenameini,status='old') snapshot=0 - DO - READ(unitini,'(a)',iostat=ios,END=100)fileheadini - IF(ios<0)EXIT ! Cycle until the last recorded state - IF(oktest) WRITE(unitterm,*)'fileheadini=',fileheadini(1:30) - READ(unitini,*,iostat=ios)it,t,ndimini,neqparini,nwini - IF(oktest) WRITE(unitterm, & + do + read(unitini,'(a)',iostat=ios,end=100)fileheadini + if(ios<0)exit ! Cycle until the last recorded state + if(oktest) write(unitterm,*)'fileheadini=',fileheadini(1:30) + read(unitini,*,iostat=ios)it,t,ndimini,neqparini,nwini + if(oktest) write(unitterm, & "('it=',i7,' t=',g10.3,' ndim=',i3,' neqpar=',i3,' nw=',i3)")& it,t,ndimini,neqparini,nwini gencoord= ndimini<0 - CALL checkNdimNeqparNw(ndimini,neqparini,nwini,neqparin,nwin) - READ(unitini,*,iostat=ios)nx - IF(oktest) WRITE(unitterm,"('nx =',3i4)")nx - CALL setixGixMix(ix^L) - READ(unitini,*,iostat=ios)(eqpar(ieqpar),ieqpar=1,neqparin),& + call checkNdimNeqparNw(ndimini,neqparini,nwini,neqparin,nwin) + read(unitini,*,iostat=ios)nx + if(oktest) write(unitterm,"('nx =',3i4)")nx + call setixGixMix(ix^L) + read(unitini,*,iostat=ios)(eqpar(ieqpar),ieqpar=1,neqparin),& (eqparextra,ieqpar=neqparin+1,neqparini) - IF(oktest) WRITE(unitterm,*)eqpar - READ(unitini,'(a)',iostat=ios)varnamesini - IF(varnames=='default')varnames=varnamesini - IF(oktest) WRITE(unitterm,*)varnames + if(oktest) write(unitterm,*)eqpar + read(unitini,'(a)',iostat=ios)varnamesini + if(varnames=='default')varnames=varnamesini + if(oktest) write(unitterm,*)varnames - {DO ix^DB=ixmin^DB,ixmax^DB;} - READ(unitini,*,iostat=ios)(x(ix^D,idim),idim=1,ndim),& + {do ix^DB=ixmin^DB,ixmax^DB;} + read(unitini,*,iostat=ios)(x(ix^D,idim),idim=1,ndim),& (w(ix^D,iw),iw=1,nwin),(wextra,iw=nwin+1,nwini) - {END DO^D&\} - IF(ios/=0)THEN - WRITE(uniterr,*)'Stop: iostat=',ios - CALL die('Error in reading file') - END IF + {end do^D&\} + call calculate_x_edges(ix^L) + if(ios/=0)then + write(uniterr,*)'Stop: iostat=',ios + call die('Error in reading file') + end if snapshot=snapshot+1 - IF(snapshot==snapshotini)EXIT - END DO + if(snapshot==snapshotini)exit + end do -100 CONTINUE +100 continue - CLOSE(unitini) + close(unitini) - IF(oktest) WRITE(*,*)'x,w:',& + if(oktest) write(*,*)'x,w:',& x(ixtest^D,idimtest),w(ixtest^D,iwtest) - IF(oktest) WRITE(*,*)'x,w:',& + if(oktest) write(*,*)'x,w:',& x(ixtest^D,idimtest),w(ixtest^D,1:nw) - RETURN -END SUBROUTINE readfileini_asc + return +end subroutine readfileini_asc !============================================================================= -SUBROUTINE readfileini_bin(w) +subroutine readfileini_bin(w) ! Reads from unitini,filenameini in binary format. ! @@ -690,147 +720,251 @@ SUBROUTINE readfileini_bin(w) ! read unless snapshotini is set. ! The compatibility of initial data with internal parameters is checked. - USE constants - USE common_varibles + use constants + use common_variables - DOUBLE PRECISION:: w(ixG^T,nw) + real(kind=8):: w(ixG^T,nw) - LOGICAL:: fileexist - INTEGER:: ios ! 0 if not EOF, -1 if EOF, >0 if error - INTEGER:: ndimini,neqparini,neqparin,nwini,nwin ! values describing input data - INTEGER:: ix^L,idim,iw,ieqpar,snapshot - DOUBLE PRECISION:: eqparextra - CHARACTER(^LENNAME) :: varnamesini + logical:: fileexist + integer:: ios ! 0 if not EOF, -1 if EOF, >0 if error + integer:: ndimini,neqparini,neqparin,nwini,nwin ! values describing input data + integer:: ix^L,idim,iw,ieqpar,snapshot + real(kind=8):: eqparextra + character(^LENNAME) :: varnamesini !----------------------------------------------------------------------------- - oktest=INDEX(teststr,'readfileini')>=1 + oktest=index(teststr,'readfileini')>=1 - IF(oktest) WRITE(unitterm,*)'ReadFileIni' + if(oktest) write(unitterm,*)'ReadFileIni' - INQUIRE(file=filenameini,exist=fileexist) - IF(.NOT.fileexist) CALL die('Stop: file does not exist, filenameini='//& + inquire(file=filenameini,exist=fileexist) + if(.not.fileexist) call die('Stop: file does not exist, filenameini='//& filenameini) - OPEN(unitini,file=filenameini,status='old',form='unformatted') + open(unitini,file=filenameini,status='old',form='unformatted') snapshot=0 - DO - READ(unitini,iostat=ios) fileheadini !END=100 - IF(ios<0)EXIT ! Cycle until the last recorded state - IF(oktest) WRITE(unitterm,*)'fileheadini=',fileheadini(1:30) - READ(unitini,iostat=ios)it,t,ndimini,neqparini,nwini - IF(oktest) WRITE(unitterm, & + do + ! Read filehead + read(unitini,iostat=ios) fileheadini !END=100 + + if(ios<0)exit ! Cycle until the last recorded state + if(oktest) write(unitterm,*)'fileheadini=',fileheadini(1:30) + + ! Read params + read(unitini,iostat=ios)it,t,ndimini,neqparini,nwini + if(oktest) write(unitterm, & "('it=',i7,' t=',g10.3,' ndim=',i3,' neqpar=',i3,' nw=',i3)")& it,t,ndimini,neqparini,nwini gencoord= ndimini<0 - CALL checkNdimNeqparNw(ndimini,neqparini,nwini,neqparin,nwin) - READ(unitini,iostat=ios)nx - IF(oktest) WRITE(unitterm,"('nx =',3i4)")nx - CALL setixGixMix(ix^L) - READ(unitini,iostat=ios)(eqpar(ieqpar),ieqpar=1,neqparin),& + ! Validate parameters? + call checkNdimNeqparNw(ndimini,neqparini,nwini,neqparin,nwin) + + ! Read nx + read(unitini,iostat=ios)nx + if(oktest) write(unitterm,"('nx =',3i4)")nx + ! This set's up the global indicies based on nx and also + ! deals with the MPI indicies etc. + call setixGixMix(ix^L) + + print*, "____________________________________" + print*, "nx", nx + print*, "ix", ix^L + + ! Read eqpar + read(unitini,iostat=ios)(eqpar(ieqpar),ieqpar=1,neqparin),& (eqparextra,ieqpar=neqparin+1,neqparini) - IF(oktest) WRITE(unitterm,*)eqpar - READ(unitini,iostat=ios)varnamesini - IF(varnames=='default')varnames=varnamesini - IF(oktest) WRITE(unitterm,*)varnames + if(oktest) write(unitterm,*)eqpar + + ! Read varnamesini + read(unitini,iostat=ios)varnamesini + if(varnames=='default')varnames=varnamesini + if(oktest) write(unitterm,*)varnames + + ! Read x array + read(unitini,iostat=ios)(x(ix^S,idim),idim=1,ndim) + call calculate_x_edges(ix^L) - READ(unitini,iostat=ios)(x(ix^S,idim),idim=1,ndim) + ! Read w array ! To conform savefileout_bin we use loop for iw - DO iw=1,nwin - READ(unitini,iostat=ios)w(ix^S,iw) - END DO - IF(ios/=0)THEN - WRITE(uniterr,*)'Error in ReadFileIni: iostat=',ios - CALL die('Error in reading file') - END IF + do iw=1,nwin + read(unitini,iostat=ios)w(ix^S,iw) + end do + if(ios/=0)then + write(uniterr,*)'Error in ReadFileIni: iostat=',ios + call die('Error in reading file') + end if snapshot=snapshot+1 - IF(snapshot==snapshotini)EXIT - END DO + if(snapshot==snapshotini)exit + end do !100 CONTINUE - CLOSE(unitini) + close(unitini) - IF(oktest) WRITE(*,*)'x,w:',& + if(oktest) write(*,*)'x,w:',& x(ixtest^D,idimtest),w(ixtest^D,iwtest) - IF(oktest) WRITE(*,*)'x,w:',& + if(oktest) write(*,*)'x,w:',& x(ixtest^D,idimtest),w(ixtest^D,1:nw) - RETURN -END SUBROUTINE readfileini_bin + return +end subroutine readfileini_bin + +subroutine readfileini_gdf(w) + + use gdf, only: gdf_parameters_T, gdf_root_datasets_T, gdf_field_type_T + use hdf5, only: h5open_f, h5gopen_f, h5fopen_f, h5fclose_f, h5close_f, HID_T, H5F_ACC_RDONLY_F + use sacgdf, only: sacgdf_read_file, build_x_array, sacgdf_read_datasets + use common_variables, only: ixGlo^D, ixGhi^D, nw, filenameini, nx, x, t, gencoord, fileheadini, unitterm, teststr, x_left_edge, x_right_edge + + implicit none + + real(kind=8), intent(inout) :: w(ixG^T,nw) + + logical :: oktest + integer(kind=4) :: error + integer(HID_T) :: file_id, grid_g_id, grid_z_id, plist_id + type(gdf_parameters_T) :: gdf_sp + type(gdf_root_datasets_T) :: gdf_rd + type(gdf_field_type_T), dimension(:), allocatable :: field_types + character(len=60) :: software_name, software_version + !class(*), dimension(:, :, :), pointer :: r_ptr + real(kind=8), dimension(:, :, :), allocatable, target :: wdata3D + integer(kind=4), dimension(^ND) :: disk_nx + + integer:: ndimini,neqparini,neqparin,nwini,nwin ! values describing input data + integer:: ix^L + + oktest=index(teststr,'readfileini')>=1 + + ! just in case you are reading in a gdf and saving out a binary + fileheadini = 'gdf' + + ! Open the interface + call h5open_f(error) + + ! Open a file for reading only + if(oktest) write(unitterm,*) "Reading GDF file: ", trim(filenameini) + call h5fopen_f(trim(filenameini), H5F_ACC_RDONLY_F, file_id, error) + + ! init the objects + call gdf_sp%init() + call gdf_rd%init(1) + + ! Read the file + call sacgdf_read_file(file_id, software_name, software_version, & + gdf_rd, gdf_sp, field_types) + + ! Decode the gdf data structures + t = gdf_sp%current_time(1) + ndimini = gdf_sp%dimensionality(1) + nwini = size(field_types, 1) + ! Set neqpar + {^IFONED neqparini = 5} + {^IFTWOD neqparini = 6} + {^IFTHREED neqparini = 7} + gencoord = ndimini<0 + ! Validate parameters? + call checkNdimNeqparNw(ndimini,neqparini,nwini,neqparin,nwin) + + nx = gdf_sp%domain_dimensions(:ndimini) + disk_nx = gdf_sp%domain_dimensions(:ndimini) + t = gdf_sp%current_time(1) + + ! This set's up the global indicies based on nx and also + ! deals with the MPI indicies etc. + call setixGixMix(ix^L) + + ! Build the x array + call build_x_array(ix^L, disk_nx, gdf_sp%domain_left_edge(:ndimini), gdf_sp%domain_right_edge(:ndimini), x) + ! Save the global left and right corners + x_left_edge = gdf_sp%domain_left_edge + x_right_edge = gdf_sp%domain_right_edge + + ! Reconstruct the w array + ! Create field groups + call h5gopen_f(file_id, "data", grid_g_id, error) !Create /data + call h5gopen_f(grid_g_id, "grid_0000000000", grid_z_id, error) !Create the top grid + + call sacgdf_read_datasets(grid_z_id, plist_id, w, ix^L) + + ! Close the file and interface + call h5fclose_f(file_id, error) + call h5close_f(error) +end subroutine readfileini_gdf !============================================================================= -SUBROUTINE checkNdimNeqparNw(ndimini,neqparini,nwini,neqparin,nwin) +subroutine checkNdimNeqparNw(ndimini,neqparini,nwini,neqparin,nwin) - USE constants - USE common_varibles + use constants + use common_variables - INTEGER:: ndimini,neqparini,nwini,neqparin,nwin + integer:: ndimini,neqparini,nwini,neqparin,nwin !----------------------------------------------------------------------------- - IF(ndim/=ABS(ndimini))THEN - WRITE(*,*)'Error in ReadFileini: ndimini=',ndimini - CALL die('Incompatible dimensionalities') - ENDIF + if(ndim/=abs(ndimini))then + write(*,*)'Error in ReadFileini: ndimini=',ndimini + call die('Incompatible dimensionalities') + endif - IF(neqpar+nspecialpar/=neqparini)WRITE(*,"(a,i3,a,i3)")& + if(neqpar+nspecialpar/=neqparini)write(*,"(a,i3,a,i3)")& 'Warning in ReadFileini: number of eq.params=',neqpar,& ' /= neqparini=',neqparini - IF(nw/=nwini)WRITE(*,"(a,i3,a,i3)")& + if(nw/=nwini)write(*,"(a,i3,a,i3)")& 'Warning in ReadFileini: number of variables nw=',nw,& ' /= nwini=',nwini - IF((neqpar+nspecialpar/=neqparini.OR.nw/=nwini).AND.varnames=='default')& - CALL die('Define varnames (in &filelist for VAC, in 3rd line for VACINI)!') + if((neqpar+nspecialpar/=neqparini.or.nw/=nwini).and.varnames=='default')& + call die('Define varnames (in &filelist for VAC, in 3rd line for VACINI)!') ! The number of equation parameters and variables to be read - neqparin=MIN(neqparini,neqpar+nspecialpar) - nwin=MIN(nwini,nw) + neqparin=min(neqparini,neqpar+nspecialpar) + nwin=min(nwini,nw) - RETURN -END SUBROUTINE checkNdimNeqparNw + return +end subroutine checkNdimNeqparNw !============================================================================= -SUBROUTINE setixGixMix(ix^L) +subroutine setixGixMix(ix^L) - USE constants - USE common_varibles + use constants + use common_variables - INTEGER:: ix^L,qnx^IFMPI + integer:: ix^L,qnx^IFMPI !----------------------------------------------------------------------------- ixGmin^D=ixGlo^D; ixMmin^D=ixGmin^D+dixBmin^D; ! Shave off ghost cells from nx - IF(fullgridini)THEN + if(fullgridini)then {^DLOOP - IF(ipe^D==0)^IFMPI nx(^D)=nx(^D)-dixBmin^D - IF(ipe^D==npe^D-1)^IFMPI nx(^D)=nx(^D)-dixBmax^D + if(ipe^D==0)^IFMPI nx(^D)=nx(^D)-dixBmin^D + if(ipe^D==npe^D-1)^IFMPI nx(^D)=nx(^D)-dixBmax^D \} - ENDIF + endif ! Calculate mesh and grid sizes ixMmax^D=ixMmin^D+nx(^D)-1; ixGmax^D=ixMmax^D+dixBmax^D; ! Set the index range for this grid - IF(fullgridini)THEN + if(fullgridini)then ix^L=ixG^L; {^IFMPI ! Set index range to mesh value if the boundary is not an outer bounary - ^D&IF(ipe^D>0)ixmin^D=ixMmin^D\ - ^D&IF(ipe^D0)ixmin^D=ixMmin^D\ + ^D&if(ipe^DixGhi^D|.OR.)THEN - WRITE(uniterr,*)'Stop: nxhi=',ixGhi^D-dixBmax^D-ixMmin^D+1 - CALL die('Error in SetixGixMix') - END IF + if(ixGmax^D>ixGhi^D|.or.)then + write(uniterr,*)'Stop: nxhi=',ixGhi^D-dixBmax^D-ixMmin^D+1 + call die('Error in SetixGixMix') + end if nx^D=nx(^D); @@ -838,225 +972,227 @@ SUBROUTINE setixGixMix(ix^L) ! set global grid size by adding up nx and ! dividing by the number of processors in the orthogonal plane {^DLOOP - CALL MPI_allreduce(nx^D,nxall^D,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierrmpi) + call MPI_allreduce(nx^D,nxall^D,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierrmpi) nxall^D=nxall^D*npe^D/npe; \} ! Distribute global grid onto processor cube - CALL mpigridsetup + call mpigridsetup } - RETURN -END SUBROUTINE setixGixMix + return +end subroutine setixGixMix !============================================================================= -SUBROUTINE setheaderstrings +subroutine setheaderstrings ! Check and/or put physics and equation parameter names into file header - USE constants - USE common_varibles + use constants + use common_variables - INTEGER:: i - CHARACTER(^LENTYPE) :: physics + integer:: i + character(^LENTYPE) :: physics !----------------------------------------------------------------------------- ! Check physics or add _typephysNDIMNDIR - WRITE(physics,'(a,i1,i1)')typephys,^ND,^NC - - i=INDEX(fileheadout,'_') - IF(i>=1)THEN - IF(physics/=fileheadout(i+1:i+^LENTYPE))THEN - WRITE(*,*)'This code is configured to ',physics - CALL die('Error: physics in file is '//fileheadout(i+1:i+^LENTYPE)) - ENDIF - ELSE - i=^LENNAME-LEN(typephys)-3 - DO - IF(fileheadout(i:i)/=' ' .OR. i==1)EXIT + write(physics,'(a,i1,i1)')typephys,^ND,^NC + + i=index(fileheadout,'_') + if(i>=1)then + if(physics/=fileheadout(i+1:i+^LENTYPE))then + write(*,*)'This code is configured to ',physics + call die('Error: physics in file is '//fileheadout(i+1:i+^LENTYPE)) + endif + else + i=^LENNAME-len(typephys)-3 + do + if(fileheadout(i:i)/=' ' .or. i==1)exit i=i-1 - ENDDO + enddo fileheadout=fileheadout(1:i)//'_'//physics - WRITE(*,*)'Warning: incomplete input headline.',& + write(*,*)'Warning: incomplete input headline.',& ' Added to output headline _',physics - ENDIF + endif ! Check for equation parameter names in varnames, add them if missing - IF(varnames/='default' .AND. INDEX(varnames,eqparname)<=0)THEN - i=^LENNAME-LEN(eqparname)-3 - DO - IF(varnames(i:i)/=' ' .OR. i==1)EXIT + if(varnames/='default' .and. index(varnames,eqparname)<=0)then + i=^LENNAME-len(eqparname)-3 + do + if(varnames(i:i)/=' ' .or. i==1)exit i=i-1 - ENDDO + enddo varnames=varnames(1:i)//' '//eqparname - ENDIF + endif ! Check for special equation parameter names in varnames, add them if missing - IF(varnames/='default' .AND. INDEX(varnames,specialparname)<=0)THEN - i=^LENNAME-LEN(specialparname)-3 - DO - IF(varnames(i:i)/=' ' .OR. i==1)EXIT + if(varnames/='default' .and. index(varnames,specialparname)<=0)then + i=^LENNAME-len(specialparname)-3 + do + if(varnames(i:i)/=' ' .or. i==1)exit i=i-1 - ENDDO + enddo varnames=varnames(1:i)//' '//specialparname - ENDIF + endif - RETURN -END SUBROUTINE setheaderstrings + return +end subroutine setheaderstrings !============================================================================= -SUBROUTINE savefile(ifile,w) +subroutine savefile(ifile,w) - USE constants - USE common_varibles + use constants + use common_variables - INTEGER:: ifile,ix^L - DOUBLE PRECISION:: w(ixG^T,nw) - CHARACTER(10):: itstring + integer:: ifile,ix^L + real(kind=8):: w(ixG^T,nw) + character(10):: itstring !----------------------------------------------------------------------------- - + print*, "in savefile" !if(nproc(ifile+2)>0) call process(ifile+2,1,ndim,w) ! In most cases the mesh should be saved ix^L=ixM^L; - SELECT CASE(ifile) - CASE(fileout_) + select case(ifile) + case(fileout_) ! Produce the output file name filenameout=filename(fileout_) - IF(snapshotout>0.AND.isaveout>0)THEN - IF(isaveout==snapshotout*(isaveout/snapshotout))THEN - CLOSE(unitini+ifile) - WRITE(itstring,'(i10)')isaveout+1 - filenameout=filenameout(1:INDEX(filenameout,' ')-1)//'_'// & - itstring(10-INT(alog10(isaveout+1.5)):10) - ENDIF - ENDIF + if(snapshotout>0.and.isaveout>0)then + if(isaveout==snapshotout*(isaveout/snapshotout))then + close(unitini+ifile) + write(itstring,'(i10)')isaveout+1 + filenameout=filenameout(1:index(filenameout,' ')-1)//'_'// & + itstring(10-int(alog10(isaveout+1.5)):10) + endif + endif isaveout=isaveout+1 - IF(fullgridout)THEN + if(fullgridout)then ix^L=ixG^L; {^IFMPI ! Set index range to mesh value if not an outer boundary - ^D&IF(ipe^D>0)ixmin^D=ixMmin^D\ - ^D&IF(ipe^D0)ixmin^D=ixMmin^D\ + ^D&if(ipe^D5.0d-16) + where(abs(w(ix^D,1:nw))>5.0d-16) qw(1:nw)=w(ix^D,1:nw) - ELSEWHERE + elsewhere qw(1:nw)=0d0 - END WHERE - WRITE(qunit,"(100(1pe18.10))")x(ix^D,1:ndim),qw(1:nw) -ENDDO^D&; + end where + write(qunit,"(100(1pe18.10))")x(ix^D,1:ndim),qw(1:nw) +enddo^D&; -CALL flushunit(qunit) +call flushunit(qunit) -RETURN -END SUBROUTINE savefileout_asc +return +end subroutine savefileout_asc !============================================================================= -SUBROUTINE savefileout_bin(qunit,w,ix^L) +subroutine savefileout_bin(qunit,w,ix^L) ! This version saves into filename(fileout_) binary data at every save time ! in full accordance with the ReadFileini subroutine, except that the first ! line is fileheadout and not fileheadini. - USE constants - USE common_varibles + use constants + use common_variables - INTEGER:: qunit,ix^L,idim,iw,ndimout - DOUBLE PRECISION:: w(ixG^T,nw) - LOGICAL:: fileopen + integer:: qunit,ix^L,idim,iw,ndimout + real(kind=8):: w(ixG^T,nw) + logical:: fileopen !**************** slice - INTEGER:: s_ixmax^D, prom_ixmax^D, s_ixmin^D, prom_ixmin^D + integer:: s_ixmax^D, prom_ixmax^D, s_ixmin^D, prom_ixmin^D !**************** endslice !----------------------------------------------------------------------------- - INQUIRE(qunit,opened=fileopen) - IF(.NOT.fileopen)& - OPEN(qunit,file=filenameout,status='unknown',form='unformatted') + inquire(qunit,opened=fileopen) + if(.not.fileopen)& + open(qunit,file=filenameout,status='unknown',form='unformatted') - IF(gencoord)THEN + if(gencoord)then ndimout= -ndim - ELSE + else ndimout= ndim - ENDIF + endif - WRITE(qunit)fileheadout - WRITE(qunit)it,t,ndimout,neqpar+nspecialpar,nw - WRITE(qunit) ixmax^D-ixmin^D+1 - WRITE(qunit)eqpar - WRITE(qunit)varnames - WRITE(qunit)(x(ix^S,idim),idim=1,ndim) + write(qunit)fileheadout + write(qunit)it,t,ndimout,neqpar+nspecialpar,nw + write(qunit) ixmax^D-ixmin^D+1 + write(qunit)eqpar + write(qunit)varnames + write(qunit)(x(ix^S,idim),idim=1,ndim) ! write(qunit)w(ix^S,1:nw) produces segmentation fault on Alpha, thus loop - DO iw=1,nw - WRITE(qunit)w(ix^S,iw) - END DO + do iw=1,nw + write(qunit)w(ix^S,iw) + end do - CALL flushunit(qunit) + call flushunit(qunit) !**************** slice ********************************* @@ -1094,11 +1230,120 @@ SUBROUTINE savefileout_bin(qunit,w,ix^L) !************* end slice ********************************** - RETURN -END SUBROUTINE savefileout_bin + return +end subroutine savefileout_bin + +!============================================================================= + + +subroutine savefileout_gdf(w,ix^L) + + use hdf5 + use gdf + use sacgdf + use gdf_datasets, only: write_dataset + use common_variables + + implicit none + + integer, intent(IN) :: ix^L + real(kind=8), intent(IN):: w(ixG^T,nw) + + integer(HID_T) :: file_id + integer(HID_T) :: plist_id, xfer_prp !< Property list identifier + integer(HID_T) :: doml_g_id !< domain list identifier + integer(HID_T) :: dom_g_id !< domain group identifier + integer :: error + + character(len=8) :: itstr + + integer, dimension(3) :: gdf_nx + type(gdf_root_datasets_T) :: rd + type(gdf_parameters_T) :: gdf_sp + type(gdf_field_type_T), dimension(nw) :: field_types + + class(*), dimension(:, :, :), pointer :: d_ptr + real(kind=8), dimension(ix^S) :: wdata + real(kind=8), dimension(:, :, :), allocatable, target :: wdata3D + real(kind=8), dimension(ixmax1, ixmax2, 2) :: temp_x + integer(kind=8), dimension(3) :: offset, count, file_dims + + gdf_nx = (/ 1, 1, 1 /) + gdf_nx(:^ND) = (/ ixmax^D-ixmin^D+1 /) + file_dims = gdf_nx + {^IFMPI {file_dims(^D) = (ixmax^D-ixmin^D+1) * npe^D; }} + + allocate(wdata3D(1:gdf_nx(1), 1:gdf_nx(2), 1:gdf_nx(3))) + + ! Open file + call h5open_f(error) + ! Create a property access list (for MPI later on) + call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, error) + ! Create the file + ! Convert the current iteration to a string + write(itstr, '(I8.8)') int(it) + call h5fcreate_f(trim(filenameout)//itstr//'.gdf', H5F_ACC_TRUNC_F, file_id, error, access_prp=plist_id) + + ! Create + ! Simulation Parameters + call gdf_sp%init() + gdf_sp%boundary_conditions = (/ 2, 2, 2, 2, 2, 2 /) + gdf_sp%cosmological_simulation = 0 + gdf_sp%current_time = t + gdf_sp%dimensionality = ^ND + gdf_sp%domain_dimensions = file_dims ! on disk + gdf_sp%domain_left_edge = x_left_edge + gdf_sp%domain_right_edge = x_right_edge + gdf_sp%field_ordering = 1 + gdf_sp%num_ghost_zones = 0 !on disk + gdf_sp%refine_by = 0 + gdf_sp%unique_identifier = "sacgdf2014" + + ! Initilize the data + call rd%init(1) + + rd%grid_parent_id = 0 + rd%grid_left_index(:, 1) = (/ 0, 0 /) + rd%grid_dimensions(:, 1) = file_dims + rd%grid_level = 0 + rd%grid_particle_count(:, 1) = (/ 0 /) + + call sacgdf_make_field_types(field_types) + + call sacgdf_write_file(file_id, rd, gdf_sp, field_types) + + ! Now write the datasets + ! Create field groups + call h5gcreate_f(file_id, "data", dom_g_id, error) !Create /data + call h5gcreate_f(dom_g_id, "grid_0000000000", doml_g_id, error) !Create the top grid + + ! WRITE ACTUAL DATA HERE + + !Calculate offset and count + offset = (/ 0, 0, 0 /) + count = (/ 1, 1, 1 /) + {count(^D) = ixmax^D-ixmin^D+1; } + ! If we are not in MPI mode use the default xfer_prp + call h5pcreate_f(H5P_DATASET_XFER_F, xfer_prp, error) + + {^IFMPI + {offset(^D) = (ixmax^D-ixmin^D+1) * ipe^D; } + {count(^D) = ixmax^D-ixmin^D+1; } + + ! Create property list for collective dataset write + !call h5pcreate_f(H5P_DATASET_XFER_F, xfer_prp, error) + !call h5pset_dxpl_mpio_f(xfer_prp, H5FD_MPIO_COLLECTIVE_F, error) + } + + call sacgdf_write_datasets(doml_g_id, w, ix^L, xfer_prp, file_dims, offset, count) + + call h5fclose_f(file_id, error) + call h5close_f(error) + +end subroutine savefileout_gdf !============================================================================= -SUBROUTINE savefilelog_default(qunit,w,ix^L) +subroutine savefilelog_default(qunit,w,ix^L) ! This version saves into filename(filelog_) the following formatted data: ! @@ -1112,67 +1357,67 @@ SUBROUTINE savefilelog_default(qunit,w,ix^L) ! at every save time. wmean is the volume averaged w, residual is saved ! if residmin>0 is set in the parfile. - USE constants - USE common_varibles + use constants + use common_variables - INTEGER:: qunit,ix^L - DOUBLE PRECISION:: w(ixG^T,nw) - INTEGER:: iw - LOGICAL:: fileopen - DOUBLE PRECISION:: wmean(nw) + integer:: qunit,ix^L + real(kind=8):: w(ixG^T,nw) + integer:: iw + logical:: fileopen + real(kind=8):: wmean(nw) !----------------------------------------------------------------------------- - IF(ipe==0)THEN^IFMPI + if(ipe==0)then^IFMPI - INQUIRE(qunit,opened=fileopen) - IF(.NOT.fileopen)THEN - OPEN(qunit,file=filename(filelog_),status='unknown') - WRITE(qunit,'(a)')fileheadout - IF(residmin>zero.OR.residmaxzero.or.residmaxzero.OR.residmaxzero.or.residmax wdata3D + call write_dataset(place, 'velocity_x', d_ptr, xfer_prp, file_dims, count, offset) + ! Density pert + wdata(ix^S) = w(ix^S, rho_) + wdata3D = reshape(wdata, gdf_nx) + d_ptr => wdata3D + call write_dataset(place, 'density_pert', d_ptr, xfer_prp, file_dims, count, offset) + + ! Denisty bg + wdata(ix^S) = w(ix^S, rhob_) + wdata3D = reshape(wdata, gdf_nx) + d_ptr => wdata3D + call write_dataset(place, 'density_bg', d_ptr, xfer_prp, file_dims, count, offset) + + ! Mag field pert + wdata(ix^S) = w(ix^S, b1_)*sqrt(mu0) + wdata3D = reshape(wdata, gdf_nx) + d_ptr => wdata3D + call write_dataset(place, 'mag_field_x_pert', d_ptr, xfer_prp, file_dims, count, offset) + + ! Mag field bg + wdata(ix^S) = w(ix^S, bg1_)*sqrt(mu0) + wdata3D = reshape(wdata, gdf_nx) + d_ptr => wdata3D + call write_dataset(place, 'mag_field_x_bg', d_ptr, xfer_prp, file_dims, count, offset) + + ! internal energy pert + wdata(ix^S) = w(ix^S, e_) + wdata3D = reshape(wdata, gdf_nx) + d_ptr => wdata3D + call write_dataset(place, 'internal_energy_pert', d_ptr, xfer_prp, file_dims, count, offset) + + ! internal energy bg + wdata(ix^S) = w(ix^S, eb_) + wdata3D = reshape(wdata, gdf_nx) + d_ptr => wdata3D + call write_dataset(place, 'internal_energy_bg', d_ptr, xfer_prp, file_dims, count, offset) + + + end subroutine sacgdf_write_datasets_1D + + +{^IFTWOD + subroutine sacgdf_write_datasets_2D(place, w, ix^L, xfer_prp, file_dims, offset, count) + use hdf5, only: HID_T + use gdf_datasets, only: write_dataset + use common_variables, only: rho_, rhob_, m2_, b2_, bg2_ + use common_variables, only: ixGhi^D, ixGlo^D, nw, nx + use phys_constants, only: mu0 + + implicit none + + integer(HID_T), intent(inout) :: place + integer(HID_T), intent(inout) :: xfer_prp + integer(kind=8), dimension(3), intent(inout) :: offset, count, file_dims + integer, intent(in) :: ix^L + real(kind=8), intent(in) :: w(ixG^T,nw) + + integer :: error + integer, dimension(3) :: gdf_nx + class(*), dimension(:, :, :), pointer :: d_ptr + real(kind=8), dimension(ix^S), target :: wdata + real(kind=8), dimension(:, :, :), allocatable, target :: wdata3D + + gdf_nx = (/ 1, 1, 1 /) + gdf_nx(:^ND) = (/ ixmax^D-ixmin^D+1 /) + allocate(wdata3D(1:gdf_nx(1), 1:gdf_nx(2), 1:gdf_nx(3))) + + ! Velocity + wdata(ix^S) = w(ix^S, m2_) / (w(ix^S, rho_) + w(ix^S, rhob_)) + wdata3D = reshape(wdata, gdf_nx) + d_ptr => wdata3D + call write_dataset(place, 'velocity_y', d_ptr, xfer_prp, file_dims, count, offset) + + ! Mag field pert + wdata(ix^S) = w(ix^S, b2_)*sqrt(mu0) + wdata3D = reshape(wdata, gdf_nx) + d_ptr => wdata3D + call write_dataset(place, 'mag_field_y_pert', d_ptr, xfer_prp, file_dims, count, offset) + + ! Mag field bg + wdata(ix^S) = w(ix^S, bg2_)*sqrt(mu0) + wdata3D = reshape(wdata, gdf_nx) + d_ptr => wdata3D + call write_dataset(place, 'mag_field_y_bg', d_ptr, xfer_prp, file_dims, count, offset) + + + end subroutine sacgdf_write_datasets_2D +} +{^IFTHREED + subroutine sacgdf_write_datasets_3D(place, w, ix^L, xfer_prp, file_dims, offset, count) + use hdf5, only: HID_T + use gdf_datasets, only: write_dataset + use common_variables, only: rho_, rhob_, m3_, b3_, bg3_ + use common_variables, only: ixGhi^D, ixGlo^D, nw, nx + use phys_constants, only: mu0 + + implicit none + + integer(HID_T), intent(inout) :: place + integer(HID_T), intent(inout) :: xfer_prp + integer(kind=8), dimension(3), intent(inout) :: offset, count, filedims + integer, intent(in) :: ix^L + real(kind=8), intent(in) :: w(ixG^T,nw) + + integer :: error + integer, dimension(3) :: gdf_nx + class(*), dimension(:, :, :), pointer :: d_ptr + real(kind=8), dimension(ix^S), target :: wdata + real(kind=8), dimension(:, :, :), allocatable, target :: wdata3D + + gdf_nx = (/ 1, 1, 1 /) + gdf_nx(:^ND) = (/ ixmax^D-ixmin^D+1 /) + allocate(wdata3D(1:gdf_nx(1), 1:gdf_nx(2), 1:gdf_nx(3))) + + ! Velocity + wdata(ix^S) = w(ix^S, m2_) / (w(ix^S, rho_) + w(ix^S, rhob_)) + wdata3D = reshape(wdata, gdf_nx) + d_ptr => wdata3D + call write_dataset(place, 'velocity_y', d_ptr, xfer_prp, file_dims, count, offset) + + ! Mag field pert + wdata(ix^S) = w(ix^S, b2_)*sqrt(mu0) + wdata3D = reshape(wdata, gdf_nx) + d_ptr => wdata3D + call write_dataset(place, 'mag_field_y_pert', d_ptr, xfer_prp, file_dims, count, offset) + + ! Mag field bg + wdata(ix^S) = w(ix^S, bg2_)*sqrt(mu0) + wdata3D = reshape(wdata, gdf_nx) + d_ptr => wdata3D + call write_dataset(place, 'mag_field_y_bg', d_ptr, xfer_prp, file_dims, count, offset) + + ! Velocity + wdata(ix^S) = w(ix^S, m3_) / (w(ix^S, rho_) + w(ix^S, rhob_)) + wdata3D = reshape(wdata, gdf_nx) + d_ptr => wdata3D + call write_dataset(place, 'velocity_z', d_ptr, xfer_prp, file_dims, count, offset) + + ! Mag field pert + wdata(ix^S) = w(ix^S, b3_)*sqrt(mu0) + wdata3D = reshape(wdata, gdf_nx) + d_ptr => wdata3D + call write_dataset(place, 'mag_field_z_pert', d_ptr, xfer_prp, file_dims, count, offset) + + ! Mag field bg + wdata(ix^S) = w(ix^S, bg3_)*sqrt(mu0) + wdata3D = reshape(wdata, gdf_nx) + d_ptr => wdata3D + call write_dataset(place, 'mag_field_z_bg', d_ptr, xfer_prp, file_dims, count, offset) + + + end subroutine sacgdf_write_datasets_3D +} + +!================================================================================ +!================================================================================ diff --git a/sac/src/vacmpi.t b/sac/src/vacmpi.t index 3ddc131..478c75f 100644 --- a/sac/src/vacmpi.t +++ b/sac/src/vacmpi.t @@ -1,177 +1,187 @@ !============================================================================= -SUBROUTINE mpiinit +subroutine mpiinit ! Initialize MPI variables - USE constants - USE common_varibles + use constants + use common_variables !---------------------------------------------------------------------------- - CALL MPI_INIT(ierrmpi) - CALL MPI_COMM_RANK (MPI_COMM_WORLD, ipe, ierrmpi) - CALL MPI_COMM_SIZE (MPI_COMM_WORLD, npe, ierrmpi) + call MPI_INIT(ierrmpi) + call MPI_COMM_RANK (MPI_COMM_WORLD, ipe, ierrmpi) + call MPI_COMM_SIZE (MPI_COMM_WORLD, npe, ierrmpi) ! unset values for directional processor numbers npe^D=-1; ! default value for test processor ipetest=0 - RETURN -END SUBROUTINE mpiinit + return +end subroutine mpiinit !============================================================================== -SUBROUTINE mpifinalize +subroutine mpifinalize - USE constants - USE common_varibles + use constants + use common_variables - CALL MPI_BARRIER(MPI_COMM_WORLD,ierrmpi) - CALL MPI_FINALIZE(ierrmpi) + call MPI_BARRIER(MPI_COMM_WORLD,ierrmpi) + call MPI_FINALIZE(ierrmpi) - RETURN -END SUBROUTINE mpifinalize + return +end subroutine mpifinalize !============================================================================== -SUBROUTINE ipe2ipeD(qipe,qipe^D) +subroutine ipe2ipeD(qipe,qipe^D) ! Convert serial processor index to directional processor indexes - USE constants - USE common_varibles + use constants + use common_variables - INTEGER:: qipe^D, qipe + integer:: qipe^D, qipe !----------------------------------------------------------------------------- qipe1 = qipe - npe1*(qipe/npe1) {qipe2 = qipe/npe1 - npe2*(qipe/(npe1*npe2)) ^NOONED} {qipe3 = qipe/(npe1*npe2) ^IFTHREED} - RETURN -END SUBROUTINE ipe2ipeD + return +end subroutine ipe2ipeD !============================================================================== -SUBROUTINE ipeD2ipe(qipe^D,qipe) +subroutine ipeD2ipe(qipe^D,qipe) ! Convert directional processor indexes to serial processor index - USE constants - USE common_varibles + use constants + use common_variables - INTEGER:: qipe^D, qipe + integer:: qipe^D, qipe !----------------------------------------------------------------------------- qipe = qipe1 {^NOONED + npe1*qipe2} {^IFTHREED + npe1*npe2*qipe3} - RETURN -END SUBROUTINE ipeD2ipe + return +end subroutine ipeD2ipe !============================================================================== -SUBROUTINE mpisetnpeDipeD(name) +subroutine mpisetnpeDipeD(name, filetype) ! Set directional processor numbers and indexes based on a filename. ! The filename contains _np followed by np^D written with 2 digit integers. ! For example _np0203 means np1=2, np2=3 for 2D. - - USE constants - USE common_varibles - CHARACTER(^LENNAME) :: name, nametail - INTEGER:: i,qnpe^D - LOGICAL:: npeDknown,npeDinname + use constants + use common_variables + + implicit none + + character(^LENNAME), intent(inout) :: name + character(^LENTYPE), intent(in) :: filetype + character(^LENNAME) :: nametail + integer:: i,qnpe^D + logical:: npeDknown,npeDinname !----------------------------------------------------------------------------- - oktest = INDEX(teststr,'mpisetnpeDipeD')>0 - IF(oktest)WRITE(*,*)'mpisetnpeDipeD ipe,name=',ipe,name + oktest = index(teststr,'mpisetnpeDipeD')>0 + if(oktest)write(*,*)'mpisetnpeDipeD ipe,name=',ipe,name ! Check if npe^D is already known npeDknown = npe1>0 - IF(npedknown .AND. npe^D* /= npe)THEN - WRITE(*,*)'npe=',npe,' /= product of npe^D=',npe^D - CALL mpistop('ERROR in setnpeDipeD') - ENDIF + if(npedknown .and. npe^D* /= npe)then + write(*,*)'npe=',npe,' /= product of npe^D=',npe^D + call mpistop('ERROR in setnpeDipeD') + endif ! Check if npe^D is given in the name - i=INDEX(name,'_np')+3 + i=index(name,'_np')+3 npeDinname = i>3 - IF(.NOT.(npeDknown.OR.npeDinname))CALL mpistop( & + if(.not.(npeDknown.or.npeDinname))call mpistop( & 'ERROR in setnpeDipeD: npeD is neither known nor given in name='//name) - IF(npeDinname)THEN + if(npeDinname)then ! read npe^D from name - READ(name(i:i+5),'(3i2)') qnpe^D + read(name(i:i+5),'(3i2)') qnpe^D i=i+2*^ND nametail=name(i:^LENNAME) - ENDIF + endif - IF( npeDknown .AND. npeDinname )THEN + if( npeDknown .and. npeDinname )then ! Check agreement - IF( qnpe^D/=npe^D|.OR. )THEN - WRITE(*,*)'npe^D=',npe^D,' /= qnpe^D=',qnpe^D,' read from filename=',name - CALL mpistop('ERROR in mpisetnpeDipeD') - ENDIF - ENDIF + if( qnpe^D/=npe^D|.or. )then + write(*,*)'npe^D=',npe^D,' /= qnpe^D=',qnpe^D,' read from filename=',name + call mpistop('ERROR in mpisetnpeDipeD') + endif + endif - IF(npeDinname .AND. .NOT.npeDknown)THEN + if(npeDinname .and. .not.npeDknown)then ! set npe^D based on name npe^D=qnpe^D; - IF( npe^D* /= npe)THEN - WRITE(*,*)'npe=',npe,' /= product of npe^D=',npe^D,& + if( npe^D* /= npe)then + write(*,*)'npe=',npe,' /= product of npe^D=',npe^D,& ' read from filename=',name - CALL mpistop('ERROR in setnpeDipeD') - ENDIF - ENDIF + call mpistop('ERROR in setnpeDipeD') + endif + endif ! Get directional processor indexes - CALL ipe2ipeD(ipe,ipe^D) + call ipe2ipeD(ipe,ipe^D) - IF(npeDknown .AND. .NOT.npeDinname)THEN + if(npeDknown .and. .not.npeDinname)then ! insert npe^D into name - i=INDEX(name,'.') + i=index(name,'.') nametail=name(i:^LENNAME) - WRITE(name(i:^LENNAME),"('_np',3i2.2)") npe^D + write(name(i:^LENNAME),"('_np',3i2.2)") npe^D i = i+3+2*^ND - ENDIF + endif + + ! insert ipe number into the filename, only if not using gdf files. + if ((filetype .eq. 'outfile') .and. (.not. typefileout .eq. 'gdf')) then + write(name(i:^LENNAME),"('_',i3.3,a)") ipe,nametail(1:^LENNAME-i-4) + end if - ! insert ipe number into the filename - WRITE(name(i:^LENNAME),"('_',i3.3,a)") ipe,nametail(1:^LENNAME-i-4) + if ((filetype .eq. 'inifile') .and. (.not. typefileini .eq. 'gdf')) then + write(name(i:^LENNAME),"('_',i3.3,a)") ipe,nametail(1:^LENNAME-i-4) + end if ! Set logicals about MPI boundaries for this processor {^DLOOP mpiupperB(^D)=ipe^D0 \} - IF(oktest)WRITE(*,*)'mpisetnpeDipeD: ipe,npeD,ipeD,name=',ipe,npe^D,ipe^D,name + if(oktest)write(*,*)'mpisetnpeDipeD: ipe,npeD,ipeD,name=',ipe,npe^D,ipe^D,name - RETURN -END SUBROUTINE mpisetnpeDipeD + return +end subroutine mpisetnpeDipeD !============================================================================== -SUBROUTINE mpineighbors(idir,hpe,jpe) +subroutine mpineighbors(idir,hpe,jpe) ! Find the hpe and jpe processors on the left and right side of this processor ! in direction idir. The processor cube is taken to be periodic in every ! direction. - USE constants - USE common_varibles + use constants + use common_variables - INTEGER :: idir,hpe,jpe,hpe^D,jpe^D + integer :: idir,hpe,jpe,hpe^D,jpe^D !----------------------------------------------------------------------------- hpe^D=ipe^D-kr(^D,idir); jpe^D=ipe^D+kr(^D,idir); {^DLOOP - IF(hpe^D<0)hpe^D=npe^D-1 - IF(jpe^D>=npe^D)jpe^D=0\} + if(hpe^D<0)hpe^D=npe^D-1 + if(jpe^D>=npe^D)jpe^D=0\} - CALL ipeD2ipe(hpe^D,hpe) - CALL ipeD2ipe(jpe^D,jpe) + call ipeD2ipe(hpe^D,hpe) + call ipeD2ipe(jpe^D,jpe) - RETURN -END SUBROUTINE mpineighbors + return +end subroutine mpineighbors !============================================================================== -SUBROUTINE mpigridsetup +subroutine mpigridsetup ! Distribute a grid of size nxall^D onto PE-s arranged in a cube of size npe^D - USE constants - USE common_varibles + use constants + use common_variables !----------------------------------------------------------------------------- !!!write(*,*)'nxall,npe=',nxall^D,npe^D @@ -184,62 +194,62 @@ SUBROUTINE mpigridsetup ! The last processors in a direction may have smaller grid sizes than nxpe {^DLOOP - IF(ipe^D < npe^D-1)THEN + if(ipe^D < npe^D-1)then nx^D = nxpe^D - ELSE + else nx^D = nxall^D - ixpemin^D + 1 - ENDIF + endif \} ! Global grid index of the last grid point stored on this PE ixPEmax^D=ixPEmin^D+nx^D-1; - RETURN -END SUBROUTINE mpigridsetup + return +end subroutine mpigridsetup !============================================================================= -SUBROUTINE mpireduce(a,mpifunc) +subroutine mpireduce(a,mpifunc) ! reduce input for one PE 0 using mpifunc - USE constants + use constants - DOUBLE PRECISION :: a, alocal - INTEGER :: mpifunc, ierrmpi + real(kind=8) :: a, alocal + integer :: mpifunc, ierrmpi !---------------------------------------------------------------------------- alocal = a - CALL MPI_REDUCE(alocal,a,1,MPI_DOUBLE_PRECISION,mpifunc,& + call MPI_REDUCE(alocal,a,1,MPI_DOUBLE_PRECISION,mpifunc,& 0,MPI_COMM_WORLD,ierrmpi) - RETURN -END SUBROUTINE mpireduce + return +end subroutine mpireduce !============================================================================== -SUBROUTINE mpiallreduce(a,mpifunc) +subroutine mpiallreduce(a,mpifunc) ! reduce input onto all PE-s using mpifunc - USE constants + use constants - DOUBLE PRECISION :: a, alocal - INTEGER :: mpifunc, ierrmpi + real(kind=8) :: a, alocal + integer :: mpifunc, ierrmpi !----------------------------------------------------------------------------- alocal = a - CALL MPI_ALLREDUCE(alocal,a,1,MPI_DOUBLE_PRECISION,mpifunc,& + call MPI_ALLREDUCE(alocal,a,1,MPI_DOUBLE_PRECISION,mpifunc,& MPI_COMM_WORLD,ierrmpi) - RETURN -END SUBROUTINE mpiallreduce + return +end subroutine mpiallreduce !============================================================================== -SUBROUTINE mpiix(ix^D,jpe) +subroutine mpiix(ix^D,jpe) ! Convert ix^D physical cell index on the global grid to local indexes ! and set the processor number jpe to the processor that contains the cell - USE constants - USE common_varibles - INTEGER :: ix^D, jpe, jpe^D + use constants + use common_variables + integer :: ix^D, jpe, jpe^D !----------------------------------------------------------------------------- ! Directional processor indexes @@ -249,83 +259,83 @@ SUBROUTINE mpiix(ix^D,jpe) ix^D=ix^D-jpe^D*nxpe^D; ! Get MPI processor index - CALL ipeD2ipe(jpe^D,jpe) + call ipeD2ipe(jpe^D,jpe) - RETURN -END SUBROUTINE mpiix + return +end subroutine mpiix !============================================================================== -SUBROUTINE mpiixlimits(ix^L) +subroutine mpiixlimits(ix^L) ! Convert global index limits to local index limits for this PE - USE constants - USE common_varibles - INTEGER :: ix^L + use constants + use common_variables + integer :: ix^L !----------------------------------------------------------------------------- {^DLOOP - IF(ixmin^D > ixPEmax^D)THEN + if(ixmin^D > ixPEmax^D)then ixmin^D = nx^D ixmax^D = nx^D-1 - ELSEIF(ixmax^D < ixPEmin^D)THEN + elseif(ixmax^D < ixPEmin^D)then ixmax^D = 0 ixmin^D = 1 - ELSE - ixmin^D = MAX(ixmin^D,ixPEmin^D) - ixPEmin^D + 1 - ixmax^D = MIN(ixmax^D,ixPEmax^D) - ixPEmin^D + 1 - ENDIF + else + ixmin^D = max(ixmin^D,ixPEmin^D) - ixPEmin^D + 1 + ixmax^D = min(ixmax^D,ixPEmax^D) - ixPEmin^D + 1 + endif \} - RETURN -END SUBROUTINE mpiixlimits + return +end subroutine mpiixlimits !============================================================================== -SUBROUTINE mpistop(message) +subroutine mpistop(message) ! Stop MPI run in an orderly fashion - USE constants - USE common_varibles + use constants + use common_variables - CHARACTER(*) :: message - INTEGER :: nerrmpi + character(*) :: message + integer :: nerrmpi !------------------------------------------------------------------------------ - WRITE(*,*)'ERROR for processor',ipe,':' - WRITE(*,*)message - CALL MPI_abort(MPI_COMM_WORLD, nerrmpi, ierrmpi) + write(*,*)'ERROR for processor',ipe,':' + write(*,*)message + call MPI_abort(MPI_COMM_WORLD, nerrmpi, ierrmpi) - STOP -END SUBROUTINE mpistop + stop +end subroutine mpistop !============================================================================== -SUBROUTINE mpibound(nvar,var) +subroutine mpibound(nvar,var) ! Fill in ghost cells of var(ixG,nvar) from other processors - USE constants - USE common_varibles + use constants + use common_variables - INTEGER :: nvar - DOUBLE PRECISION :: var(ixG^T,nvar) + integer :: nvar + real(kind=8) :: var(ixG^T,nvar) ! processor indexes for left and right neighbors - INTEGER :: hpe,jpe + integer :: hpe,jpe ! index limits for the left and right side mesh and ghost cells - INTEGER :: ixLM^L, ixRM^L, ixLG^L, ixRG^L - LOGICAL :: periodic + integer :: ixLM^L, ixRM^L, ixLG^L, ixRG^L + logical :: periodic ! There can be at most 2 receives in any direction for each PE - INTEGER :: nmpirequest, mpirequests(2) - INTEGER :: mpistatus(MPI_STATUS_SIZE,2) - COMMON /mpirecv/ nmpirequest,mpirequests,mpistatus + integer :: nmpirequest, mpirequests(2) + integer :: mpistatus(MPI_STATUS_SIZE,2) + common /mpirecv/ nmpirequest,mpirequests,mpistatus !----------------------------------------------------------------------------- - oktest=INDEX(teststr,'mpibound')>0 - IF(oktest)WRITE(*,*)'mpibound ipe,nvar,varold=',& - ipe,nvar,var(ixtest^D,MIN(nvar,iwtest)) + oktest=index(teststr,'mpibound')>0 + if(oktest)write(*,*)'mpibound ipe,nvar,varold=',& + ipe,nvar,var(ixtest^D,min(nvar,iwtest)) {^DLOOP - IF(npe^D>1)THEN + if(npe^D>1)then nmpirequest =0 mpirequests(1:2) = MPI_REQUEST_NULL @@ -340,128 +350,128 @@ SUBROUTINE mpibound(nvar,var) ixRM^L=ixG^L; ixRMmax^D=ixMmax^D; ixRMmin^D=ixMmax^D-dixBmax^D+1; ! Obtain left and right neighbor processors for this direction - CALL mpineighbors(^D,hpe,jpe) + call mpineighbors(^D,hpe,jpe) - IF(oktest)THEN - WRITE(*,*)'mpibound ipe,idir=',ipe,^D - WRITE(*,*)'mpibound ipe,ixLG=',ipe,ixLG^L - WRITE(*,*)'mpibound ipe,ixRG=',ipe,ixRG^L - WRITE(*,*)'mpibound ipe,ixLM=',ipe,ixLM^L - WRITE(*,*)'mpibound ipe,ixRM=',ipe,ixRM^L - WRITE(*,*)'mpibound ipe,hpe,jpe=',ipe,hpe,jpe - ENDIF + if(oktest)then + write(*,*)'mpibound ipe,idir=',ipe,^D + write(*,*)'mpibound ipe,ixLG=',ipe,ixLG^L + write(*,*)'mpibound ipe,ixRG=',ipe,ixRG^L + write(*,*)'mpibound ipe,ixLM=',ipe,ixLM^L + write(*,*)'mpibound ipe,ixRM=',ipe,ixRM^L + write(*,*)'mpibound ipe,hpe,jpe=',ipe,hpe,jpe + endif ! receive right (2) boundary from left neighbor hpe - IF(mpilowerB(^D) .OR. periodic)CALL mpirecvbuffer(nvar,ixRM^L,hpe,2) + if(mpilowerB(^D) .or. periodic)call mpirecvbuffer(nvar,ixRM^L,hpe,2) ! receive left (1) boundary from right neighbor jpe - IF(mpiupperB(^D) .OR. periodic)CALL mpirecvbuffer(nvar,ixLM^L,jpe,1) + if(mpiupperB(^D) .or. periodic)call mpirecvbuffer(nvar,ixLM^L,jpe,1) ! Wait for all receives to be posted - CALL MPI_BARRIER(MPI_COMM_WORLD,ierrmpi) + call MPI_BARRIER(MPI_COMM_WORLD,ierrmpi) ! Ready send left (1) boundary to left neighbor hpe - IF(mpilowerB(^D) .OR. periodic)CALL mpisend(nvar,var,ixLM^L,hpe,1) + if(mpilowerB(^D) .or. periodic)call mpisend(nvar,var,ixLM^L,hpe,1) ! Ready send right (2) boundary to right neighbor - IF(mpiupperB(^D) .OR. periodic)CALL mpisend(nvar,var,ixRM^L,jpe,2) + if(mpiupperB(^D) .or. periodic)call mpisend(nvar,var,ixRM^L,jpe,2) ! Wait for messages to arrive - CALL MPI_WAITALL(nmpirequest,mpirequests,mpistatus,ierrmpi) + call MPI_WAITALL(nmpirequest,mpirequests,mpistatus,ierrmpi) ! Copy buffer received from right (2) physical cells into left ghost cells - IF(mpilowerB(^D) .OR. periodic)CALL mpibuffer2var(2,nvar,var,ixLG^L) + if(mpilowerB(^D) .or. periodic)call mpibuffer2var(2,nvar,var,ixLG^L) ! Copy buffer received from left (1) physical cells into right ghost cells - IF(mpiupperB(^D) .OR. periodic)CALL mpibuffer2var(1,nvar,var,ixRG^L) - ENDIF + if(mpiupperB(^D) .or. periodic)call mpibuffer2var(1,nvar,var,ixRG^L) + endif \} - IF(oktest)WRITE(*,*)'mpibound ipe,varnew=',ipe,var(ixtest^D,MIN(nvar,iwtest)) + if(oktest)write(*,*)'mpibound ipe,varnew=',ipe,var(ixtest^D,min(nvar,iwtest)) - RETURN -END SUBROUTINE mpibound + return +end subroutine mpibound !============================================================================== -SUBROUTINE mpisend(nvar,var,ix^L,qipe,iside) +subroutine mpisend(nvar,var,ix^L,qipe,iside) ! Send var(ix^L,1:nvar) to processor qipe. ! jside is 0 for min and 1 for max side of the grid for the sending PE - USE constants - USE common_varibles + use constants + use common_variables - INTEGER :: nvar - DOUBLE PRECISION :: var(ixG^T,nvar) - INTEGER :: ix^L, qipe, iside, n, ix^D, ivar + integer :: nvar + real(kind=8) :: var(ixG^T,nvar) + integer :: ix^L, qipe, iside, n, ix^D, ivar !---------------------------------------------------------------------------- - oktest = INDEX(teststr,'mpisend')>0 + oktest = index(teststr,'mpisend')>0 n=0 - DO ivar=1,nvar - {DO ix^DB=ixmin^DB,ixmax^DB;} + do ivar=1,nvar + {do ix^DB=ixmin^DB,ixmax^DB;} n=n+1 sendbuffer(n)=var(ix^D,ivar) - {ENDDO^DLOOP\} - END DO + {enddo^DLOOP\} + end do - IF(oktest)THEN - WRITE(*,*)'mpisend ipe-->qipe,iside,itag',ipe,qipe,iside,10*ipe+iside - WRITE(*,*)'mpisend ipe,ix^L,var=',ipe,ix^L,var(ixtest^D,MIN(iwtest,nvar)) - ENDIF + if(oktest)then + write(*,*)'mpisend ipe-->qipe,iside,itag',ipe,qipe,iside,10*ipe+iside + write(*,*)'mpisend ipe,ix^L,var=',ipe,ix^L,var(ixtest^D,min(iwtest,nvar)) + endif - CALL MPI_RSEND(sendbuffer(1),n,MPI_DOUBLE_PRECISION,qipe,10*ipe+iside,& + call MPI_RSEND(sendbuffer(1),n,MPI_DOUBLE_PRECISION,qipe,10*ipe+iside,& MPI_COMM_WORLD,ierrmpi) - RETURN -END SUBROUTINE mpisend + return +end subroutine mpisend !============================================================================== -SUBROUTINE mpirecvbuffer(nvar,ix^L,qipe,iside) +subroutine mpirecvbuffer(nvar,ix^L,qipe,iside) ! receive buffer for a ghost cell region of size ix^L sent from processor qipe ! and sent from side iside of the grid - USE constants - USE common_varibles + use constants + use common_variables - INTEGER:: nvar, ix^L, qipe, iside, n + integer:: nvar, ix^L, qipe, iside, n - INTEGER :: nmpirequest, mpirequests(2) - INTEGER :: mpistatus(MPI_STATUS_SIZE,2) - COMMON /mpirecv/ nmpirequest,mpirequests,mpistatus + integer :: nmpirequest, mpirequests(2) + integer :: mpistatus(MPI_STATUS_SIZE,2) + common /mpirecv/ nmpirequest,mpirequests,mpistatus !---------------------------------------------------------------------------- - oktest = INDEX(teststr,'mpirecv')>0 + oktest = index(teststr,'mpirecv')>0 n = nvar* ^D&(ixmax^D-ixmin^D+1)* - IF(oktest)WRITE(*,*)'mpirecv ipe<--qipe,iside,itag,n',& + if(oktest)write(*,*)'mpirecv ipe<--qipe,iside,itag,n',& ipe,qipe,iside,10*qipe+iside,n nmpirequest = nmpirequest + 1 - CALL MPI_IRECV(recvbuffer(1,iside),n,MPI_DOUBLE_PRECISION,qipe,10*qipe+iside,& + call MPI_IRECV(recvbuffer(1,iside),n,MPI_DOUBLE_PRECISION,qipe,10*qipe+iside,& MPI_COMM_WORLD,mpirequests(nmpirequest),ierrmpi) - RETURN -END SUBROUTINE mpirecvbuffer + return +end subroutine mpirecvbuffer !============================================================================== -SUBROUTINE mpibuffer2var(iside,nvar,var,ix^L) +subroutine mpibuffer2var(iside,nvar,var,ix^L) ! Copy mpibuffer(:,iside) into var(ix^L,1:nvar) - USE constants - USE common_varibles + use constants + use common_variables - INTEGER :: nvar - DOUBLE PRECISION:: var(ixG^T,nvar) - INTEGER:: ix^L,iside,n,ix^D,ivar + integer :: nvar + real(kind=8):: var(ixG^T,nvar) + integer:: ix^L,iside,n,ix^D,ivar !----------------------------------------------------------------------------- - oktest = INDEX(teststr,'buffer2var')>0 + oktest = index(teststr,'buffer2var')>0 n=0 - DO ivar=1,nvar - {DO ix^DB=ixmin^DB,ixmax^DB;} + do ivar=1,nvar + {do ix^DB=ixmin^DB,ixmax^DB;} n=n+1 var(ix^D,ivar)=recvbuffer(n,iside) - {ENDDO^DLOOP\} - END DO + {enddo^DLOOP\} + end do - IF(oktest)WRITE(*,*)'buffer2var: ipe,iside,ix^L,var',& - ipe,iside,ix^L,var(ixtest^D,MIN(iwtest,nvar)) + if(oktest)write(*,*)'buffer2var: ipe,iside,ix^L,var',& + ipe,iside,ix^L,var(ixtest^D,min(iwtest,nvar)) - RETURN -END SUBROUTINE mpibuffer2var + return +end subroutine mpibuffer2var diff --git a/sac/src/vacphys.mhd0.t b/sac/src/vacphys.mhd0.t index b221ece..1cc5e27 100644 --- a/sac/src/vacphys.mhd0.t +++ b/sac/src/vacphys.mhd0.t @@ -11,7 +11,7 @@ SUBROUTINE physini ! Tell VAC which variables are vectors, set default entropy coefficients USE constants - USE common_varibles + USE common_variables INTEGER:: il !----------------------------------------------------------------------------- @@ -41,13 +41,13 @@ SUBROUTINE process(count,idim^LIM,w) ! count=ifile+2 for saving results into the file indexed by ifile USE constants - USE common_varibles + USE common_variables INTEGER:: count,idim^LIM - DOUBLE PRECISION:: w(ixG^T,nw) + REAL(kind=8):: w(ixG^T,nw) LOGICAL:: oktime - DOUBLE PRECISION:: cputime,time1,timeproc + REAL(kind=8):: cputime,time1,timeproc DATA timeproc /0.D0/ ! The processing should eliminate divergence of B. @@ -88,9 +88,9 @@ SUBROUTINE getdt(w,ix^L) ! If resistivity is not zero, check diffusion time limit for dt USE constants - USE common_varibles + USE common_variables - DOUBLE PRECISION:: w(ixG^T,nw) + REAL(kind=8):: w(ixG^T,nw) INTEGER:: ix^L !----------------------------------------------------------------------------- @@ -112,10 +112,10 @@ SUBROUTINE getdivb(w,ixO^L,divb) ! Calculate div B within ixO USE constants - USE common_varibles + USE common_variables INTEGER:: ixO^L,ix^L,idim - DOUBLE PRECISION:: w(ixG^T,nw),divb(ixG^T) + REAL(kind=8):: w(ixG^T,nw),divb(ixG^T) !----------------------------------------------------------------------------- oktest=INDEX(teststr,'getdivb')>=1 @@ -155,10 +155,10 @@ SUBROUTINE getflux(w,ix^L,iw,idim,f,transport) ! Set transport=.true. if a transport flux should be added USE constants - USE common_varibles + USE common_variables INTEGER:: ix^L,iw,idim - DOUBLE PRECISION:: w(ixG^T,nw),f(ixG^T), fb(ixG^T) + REAL(kind=8):: w(ixG^T,nw),f(ixG^T), fb(ixG^T) LOGICAL:: transport !----------------------------------------------------------------------------- @@ -236,10 +236,10 @@ SUBROUTINE addsource(qdt,ixI^L,ixO^L,iws,qtC,w,qt,wnew) ! Add sources from resistivity and Powell solver USE constants - USE common_varibles + USE common_variables INTEGER:: ixI^L,ixO^L,iws(niw_) - DOUBLE PRECISION:: qdt,qtC,qt,w(ixG^T,nw),wnew(ixG^T,nw) + REAL(kind=8):: qdt,qtC,qt,w(ixG^T,nw),wnew(ixG^T,nw) !----------------------------------------------------------------------------- oktest=INDEX(teststr,'addsource')>=1 @@ -276,11 +276,11 @@ SUBROUTINE addsource_divb(qdt,ixI^L,ixO^L,iws,qtC,w,qt,wnew) ! otherwise shrink ixO USE constants - USE common_varibles + USE common_variables INTEGER:: ixI^L,ixO^L,iws(niw_),iiw,iw - DOUBLE PRECISION:: qdt,qtC,qt,w(ixG^T,nw),wnew(ixG^T,nw) - DOUBLE PRECISION:: divb(ixG^T) + REAL(kind=8):: qdt,qtC,qt,w(ixG^T,nw),wnew(ixG^T,nw) + REAL(kind=8):: divb(ixG^T) !----------------------------------------------------------------------------- ! Calculating div B involves first derivatives diff --git a/sac/src/vacphys.t.mhd b/sac/src/vacphys.t.mhd index 909aac4..3d867fd 100644 --- a/sac/src/vacphys.t.mhd +++ b/sac/src/vacphys.t.mhd @@ -10,10 +10,10 @@ SUBROUTINE keeppositive(ix^L,w) ! Keep pressure and density positive (following D.Ryu) USE constants - USE common_varibles + USE common_variables INTEGER:: ix^L - DOUBLE PRECISION:: w(ixG^T,nw) + REAL(kind=8):: w(ixG^T,nw) LOGICAL:: toosmallp !----------------------------------------------------------------------------- diff --git a/sac/src/vacphys0.t.mhd b/sac/src/vacphys0.t.mhd index 32dd0da..1c89835 100644 --- a/sac/src/vacphys0.t.mhd +++ b/sac/src/vacphys0.t.mhd @@ -7,10 +7,10 @@ SUBROUTINE conserve(ix^L,w) ! Transform primitive variables into conservative ones USE constants - USE common_varibles + USE common_variables INTEGER:: ix^L - DOUBLE PRECISION:: w(ixG^T,nw) + REAL(kind=8):: w(ixG^T,nw) !----------------------------------------------------------------------------- @@ -34,10 +34,10 @@ SUBROUTINE primitive(ix^L,w) ! Transform conservative variables into primitive ones USE constants - USE common_varibles + USE common_variables INTEGER:: ix^L - DOUBLE PRECISION:: w(ixG^T,nw) + REAL(kind=8):: w(ixG^T,nw) !----------------------------------------------------------------------------- @@ -59,10 +59,10 @@ SUBROUTINE getv(w,ix^L,idim,v) ! Calculate v_idim=m_idim/rho within ix USE constants - USE common_varibles + USE common_variables INTEGER:: ix^L,idim - DOUBLE PRECISION:: w(ixG^T,nw),v(ixG^T) + REAL(kind=8):: w(ixG^T,nw),v(ixG^T) !----------------------------------------------------------------------------- oktest=INDEX(teststr,'getv')>=1 @@ -85,12 +85,12 @@ SUBROUTINE getcmax(new_cmax,w,ix^L,idim,cmax) ! perpendicular to the magnetic field, and cs is the sound speed. USE constants - USE common_varibles + USE common_variables LOGICAL:: new_cmax INTEGER:: ix^L,idim - DOUBLE PRECISION:: w(ixG^T,nw),cmax(ixG^T) - DOUBLE PRECISION:: csound2(ixG^T),cfast2(ixG^T) + REAL(kind=8):: w(ixG^T,nw),cmax(ixG^T) + REAL(kind=8):: csound2(ixG^T),cfast2(ixG^T) SAVE csound2,cfast2 !----------------------------------------------------------------------------- oktest=INDEX(teststr,'getcmax')>=1 @@ -123,9 +123,9 @@ SUBROUTINE getcsound2prim(w,ix^L,csound2) ! csound2=gamma*p/rho USE constants - USE common_varibles + USE common_variables - DOUBLE PRECISION:: w(ixG^T,nw),csound2(ixG^T) + REAL(kind=8):: w(ixG^T,nw),csound2(ixG^T) INTEGER:: ix^L !----------------------------------------------------------------------------- @@ -145,14 +145,14 @@ SUBROUTINE getcsound2(w,ix^L,csound2) ! csound2=gamma*p/rho USE constants - USE common_varibles + USE common_variables - DOUBLE PRECISION:: w(ixG^T,nw),csound2(ixG^T) + REAL(kind=8):: w(ixG^T,nw),csound2(ixG^T) INTEGER:: ix^L !----------------------------------------------------------------------------- IF(eqpar(gamma_)<=zero)& - CALL die('Correct Getcsound2 for NONIDEAL gas in vacphys.t.mhd') + CALL die('FATAL: Incorrect Getcsound2 for NONIDEAL gas in vacphys.t.mhd') oktest=INDEX(teststr,'getcsound2')>=1 IF(oktest) WRITE(*,*)'Getcsound2' @@ -172,9 +172,9 @@ SUBROUTINE getpthermal(w,ix^L,p) USE constants - USE common_varibles + USE common_variables - DOUBLE PRECISION:: w(ixG^T,nw),p(ixG^T) + REAL(kind=8):: w(ixG^T,nw),p(ixG^T) INTEGER:: ix^L !----------------------------------------------------------------------------- @@ -193,9 +193,9 @@ END SUBROUTINE getpthermal SUBROUTINE getptotal(w,ix^L,p) USE constants - USE common_varibles + USE common_variables - DOUBLE PRECISION:: w(ixG^T,nw),p(ixG^T),gamma + REAL(kind=8):: w(ixG^T,nw),p(ixG^T),gamma INTEGER:: ix^L !----------------------------------------------------------------------------- @@ -216,9 +216,9 @@ END SUBROUTINE getptotal SUBROUTINE getptotal_bg(w,ix^L,p) USE constants - USE common_varibles + USE common_variables - DOUBLE PRECISION:: w(ixG^T,nw),p(ixG^T),gamma + REAL(kind=8):: w(ixG^T,nw),p(ixG^T),gamma INTEGER:: ix^L !----------------------------------------------------------------------------- diff --git a/sac/src/vacusr.gravity.t b/sac/src/vacusr.gravity.t index 1f7fa60..47768dc 100644 --- a/sac/src/vacusr.gravity.t +++ b/sac/src/vacusr.gravity.t @@ -27,10 +27,10 @@ SUBROUTINE addsource_grav(qdt,ixI^L,ixO^L,iws,qtC,w,qt,wnew) ! in iws. w is at time qtC, wnew is advanced from qt to qt+qdt. USE constants - USE common_varibles + USE common_variables INTEGER:: ixI^L,ixO^L,iws(niw_) - DOUBLE PRECISION:: qdt,qtC,qt,w(ixG^T,nw),wnew(ixG^T,nw) + REAL(kind=8):: qdt,qtC,qt,w(ixG^T,nw),wnew(ixG^T,nw) INTEGER:: iiw,iw,idim !!! ! For a spatially varying gravity define the common grav array !!! double precision:: grav(ixG^T,ndim) @@ -98,11 +98,11 @@ END SUBROUTINE addsource_grav SUBROUTINE getdt_grav(w,ix^L) USE constants - USE common_varibles + USE common_variables - DOUBLE PRECISION :: w(ixG^T,nw) + REAL(kind=8) :: w(ixG^T,nw) INTEGER :: ix^L,idim - DOUBLE PRECISION, SAVE :: dtgrav + REAL(kind=8), SAVE :: dtgrav !!! ! For spatially varying gravity you need a common grav array !!! double precision:: grav(ixG^T,ndim) diff --git a/sac/src/vacusr.t.stuart b/sac/src/vacusr.t.default similarity index 91% rename from sac/src/vacusr.t.stuart rename to sac/src/vacusr.t.default index 412cb9c..3ed077f 100644 --- a/sac/src/vacusr.t.stuart +++ b/sac/src/vacusr.t.default @@ -12,10 +12,10 @@ SUBROUTINE specialini(ix^L,w) ! Initialize w for VACINI, user-defined USE constants - USE common_varibles + USE common_variables INTEGER:: ix^L - DOUBLE PRECISION:: w(ixG^T,nw) + REAL(kind=8):: w(ixG^T,nw) !----------------------------------------------------------------------------- CALL die('Special initial condition is not defined') @@ -28,10 +28,10 @@ SUBROUTINE specialbound(qt,ix^L,iw,iB,w) ! Calculates the boundary values in the iB-th boundary segment, user-defined USE constants - USE common_varibles + USE common_variables INTEGER:: ix^L,iw,iB - DOUBLE PRECISION:: qt,w(ixG^T,nw) + REAL(kind=8):: qt,w(ixG^T,nw) !----------------------------------------------------------------------------- CALL die('Special boundary is not defined') @@ -59,10 +59,10 @@ SUBROUTINE specialsource(qdt,ixI^L,ixO^L,iws,qtC,wCT,qt,w) ! where "dix" is the number of extra layers needed, typically 1 or 2. USE constants - USE common_varibles + USE common_variables INTEGER:: ixI^L,ixO^L,iws(niw_) - DOUBLE PRECISION:: qdt,qtC,qt,wCT(ixG^T,nw),w(ixG^T,nw) + REAL(kind=8):: qdt,qtC,qt,wCT(ixG^T,nw),w(ixG^T,nw) !This is some hyperdiffusion stabilisaton... IF(ABS(eqpar(nu_))>smalldouble)& @@ -79,9 +79,9 @@ SUBROUTINE getdt_special(w,ix^L) ! module have already been called. USE constants - USE common_varibles + USE common_variables - DOUBLE PRECISION:: w(ixG^T,nw) + REAL(kind=8):: w(ixG^T,nw) INTEGER:: ix^L !----------------------------------------------------------------------------- @@ -103,9 +103,9 @@ SUBROUTINE specialeta(w,ix^L,idirmin) ! and/or local values (ix) of "current" only, i.e. it has to be compact. USE constants - USE common_varibles + USE common_variables - DOUBLE PRECISION:: w(ixG^T,nw) + REAL(kind=8):: w(ixG^T,nw) INTEGER:: ix^L,idirmin CALL die('specialeta is not defined') @@ -119,9 +119,9 @@ SUBROUTINE readfileini_special(w) ! Check readfileini_asc and readfileini_bin in vacio.t on what should be done. USE constants - USE common_varibles + USE common_variables - DOUBLE PRECISION:: w(ixG^T,nw) + REAL(kind=8):: w(ixG^T,nw) !----------------------------------------------------------------------------- CALL die('Special readfileini is not defined') @@ -133,10 +133,10 @@ SUBROUTINE savefileout_special(qunit,w,ix^L) ! Check savefileout_asc and savefileout_bin in vacio.t on what should be done. USE constants - USE common_varibles + USE common_variables INTEGER:: qunit,ix^L - DOUBLE PRECISION:: w(ixG^T,nw) + REAL(kind=8):: w(ixG^T,nw) !----------------------------------------------------------------------------- CALL die('Special savefileout is not defined') @@ -148,10 +148,10 @@ SUBROUTINE savefilelog_special(qunit,w,ix^L) ! Check savefilelog_default on opening the file etc. USE constants - USE common_varibles + USE common_variables INTEGER:: qunit,ix^L - DOUBLE PRECISION:: w(ixG^T,nw) + REAL(kind=8):: w(ixG^T,nw) !----------------------------------------------------------------------------- CALL die('Special savefilelog is not defined') diff --git a/sac/src/vacusr.viscosity.t b/sac/src/vacusr.viscosity.t index c2ec2d9..37c80cd 100644 --- a/sac/src/vacusr.viscosity.t +++ b/sac/src/vacusr.viscosity.t @@ -4,34 +4,34 @@ SUBROUTINE addsource_visc(qdt,ixI^L,ixO^L,iws,qtC,w,qt,wnew) ! Add viscosity source to wnew within ixO USE constants - USE common_varibles + USE common_variables INTEGER:: ixI^L,ixO^L,iws(niw_) - DOUBLE PRECISION:: qdt,qtC,qt,w(ixG^T,nw),wnew(ixG^T,nw) + REAL(kind=8):: qdt,qtC,qt,w(ixG^T,nw),wnew(ixG^T,nw) INTEGER:: ix,ix^L,idim,idir,jdir,iiw,iw !already declared in vacusr.f !double precision:: tmp2(ixG^T) - DOUBLE PRECISION:: nushk(ixG^T,ndim) + REAL(kind=8):: nushk(ixG^T,ndim) - DOUBLE PRECISION:: tmprhoL(ixG^T), tmprhoR(ixG^T), tmprhoC(ixG^T) - DOUBLE PRECISION:: tmpVL(ixG^T), tmpVR(ixG^T), tmpVC(ixG^T) - DOUBLE PRECISION:: tmpBL(ixG^T), tmpBR(ixG^T), tmpBC(ixG^T) + REAL(kind=8):: tmprhoL(ixG^T), tmprhoR(ixG^T), tmprhoC(ixG^T) + REAL(kind=8):: tmpVL(ixG^T), tmpVR(ixG^T), tmpVC(ixG^T) + REAL(kind=8):: tmpBL(ixG^T), tmpBR(ixG^T), tmpBC(ixG^T) - DOUBLE PRECISION:: tmpL(ixG^T),tmpR(ixG^T), tmpC(ixG^T) + REAL(kind=8):: tmpL(ixG^T),tmpR(ixG^T), tmpC(ixG^T) - DOUBLE PRECISION:: nuL(ixG^T),nuR(ixG^T) + REAL(kind=8):: nuL(ixG^T),nuR(ixG^T) INTEGER:: jx^L,hx^L, hxO^L - DOUBLE PRECISION:: c_ene,c_shk + REAL(kind=8):: c_ene,c_shk INTEGER:: i,j,k,l,m,ii0,ii1,t00 - DOUBLE PRECISION:: sB + REAL(kind=8):: sB !----------------------------------------------------------------------------- @@ -227,17 +227,17 @@ SUBROUTINE setnu(w,iw,idim,ix^L,nuR,nuL) ! Set the viscosity coefficient nu within ixO based on w(ixI). USE constants - USE common_varibles + USE common_variables INTEGER:: ixi^L - DOUBLE PRECISION:: w(ixG^T,nw) - DOUBLE PRECISION:: d1R(^SIDEADO),d1L(^SIDEADO) - DOUBLE PRECISION:: d3R(^SIDEADO),d3L(^SIDEADO) - DOUBLE PRECISION:: md3R(ixG^T),md3L(ixG^T) - DOUBLE PRECISION:: md1R(ixG^T),md1L(ixG^T) - DOUBLE PRECISION:: nuR(ixG^T),nuL(ixG^T) - - DOUBLE PRECISION:: c_tot, c_hyp,cmax(ixG^T), tmp_nu(ixG^T) + REAL(kind=8):: w(ixG^T,nw) + REAL(kind=8):: d1R(^SIDEADO),d1L(^SIDEADO) + REAL(kind=8):: d3R(^SIDEADO),d3L(^SIDEADO) + REAL(kind=8):: md3R(ixG^T),md3L(ixG^T) + REAL(kind=8):: md1R(ixG^T),md1L(ixG^T) + REAL(kind=8):: nuR(ixG^T),nuL(ixG^T) + + REAL(kind=8):: c_tot, c_hyp,cmax(ixG^T), tmp_nu(ixG^T) INTEGER:: ix^L,idim, iw INTEGER:: kx^L,jx^L,hx^L,gx^L,ixFF^L,jxFF^L,hxFF^L INTEGER:: ix_1,ix_2,ix_3 @@ -246,7 +246,7 @@ SUBROUTINE setnu(w,iw,idim,ix^L,nuR,nuL) LOGICAL:: new_cmax - DOUBLE PRECISION:: tmp_nuI(^SIDEADD) + REAL(kind=8):: tmp_nuI(^SIDEADD) INTEGER:: k,iwc @@ -261,10 +261,10 @@ SUBROUTINE setnu(w,iw,idim,ix^L,nuR,nuL) INTEGER:: hpe,jpe - DOUBLE PRECISION:: tgtbufferR^D(1^D%^LM) - DOUBLE PRECISION:: tgtbufferL^D(1^D%^LM) - DOUBLE PRECISION:: srcbufferR^D(1^D%^LM) - DOUBLE PRECISION:: srcbufferL^D(1^D%^LM) + REAL(kind=8):: tgtbufferR^D(1^D%^LM) + REAL(kind=8):: tgtbufferL^D(1^D%^LM) + REAL(kind=8):: srcbufferR^D(1^D%^LM) + REAL(kind=8):: srcbufferL^D(1^D%^LM) INTEGER:: n @@ -462,14 +462,14 @@ END SUBROUTINE setnu SUBROUTINE setnushk(w,ix^L,nushk) USE constants - USE common_varibles + USE common_variables !double precision:: w(ixG^T,nw),tmp2(ixG^T),nushk(ixG^T,ndim) - DOUBLE PRECISION:: w(ixG^T,nw),nushk(ixG^T,ndim) + REAL(kind=8):: w(ixG^T,nw),nushk(ixG^T,ndim) - DOUBLE PRECISION:: c_shk + REAL(kind=8):: c_shk - DOUBLE PRECISION:: tmp3(ixG^T) + REAL(kind=8):: tmp3(ixG^T) INTEGER:: ix^L,idim, iw,i @@ -516,15 +516,15 @@ SUBROUTINE getdt_visc(w,ix^L) ! Based on Hirsch volume 2, p.631, eq.23.2.17 USE constants - USE common_varibles + USE common_variables - DOUBLE PRECISION:: w(ixG^T,nw),dtdiff_visc + REAL(kind=8):: w(ixG^T,nw),dtdiff_visc INTEGER:: ix^L,idim, ix_1,ix_2 INTEGER:: aa ! For spatially varying nu you need a common nu array - DOUBLE PRECISION::tmpdt(ixG^T), nuL(ixG^T),nuR(ixG^T), nushk(ixG^T,ndim) + REAL(kind=8)::tmpdt(ixG^T), nuL(ixG^T),nuR(ixG^T), nushk(ixG^T,ndim) COMMON/visc/nuL COMMON/visc/nuR !----------------------------------------------------------------------------- @@ -550,9 +550,9 @@ END SUBROUTINE getdt_visc SUBROUTINE gradient1(q,ix^L,idim,gradq) USE constants - USE common_varibles + USE common_variables INTEGER:: ix^L,idim - DOUBLE PRECISION:: q(ixG^T),gradq(ixG^T) + REAL(kind=8):: q(ixG^T),gradq(ixG^T) INTEGER:: hx^L,kx^L INTEGER:: minx1^D,maxx1^D,k !----------------------------------------------------------------------------- @@ -597,9 +597,9 @@ END SUBROUTINE gradient1 SUBROUTINE gradient1L(q,ix^L,idim,gradq) USE constants - USE common_varibles + USE common_variables INTEGER:: ix^L,idim - DOUBLE PRECISION:: q(ixG^T),gradq(ixG^T) + REAL(kind=8):: q(ixG^T),gradq(ixG^T) INTEGER:: hx^L INTEGER:: minx1^D,maxx1^D,k !----------------------------------------------------------------------------- @@ -642,9 +642,9 @@ END SUBROUTINE gradient1L SUBROUTINE gradient1R(q,ix^L,idim,gradq) USE constants - USE common_varibles + USE common_variables INTEGER:: ix^L,idim - DOUBLE PRECISION:: q(ixG^T),gradq(ixG^T) + REAL(kind=8):: q(ixG^T),gradq(ixG^T) INTEGER:: hx^L INTEGER:: minx1^D,maxx1^D,k !----------------------------------------------------------------------------- diff --git a/sac/src/vacusrpar.t.stuart b/sac/src/vacusrpar.gravity.t similarity index 100% rename from sac/src/vacusrpar.t.stuart rename to sac/src/vacusrpar.gravity.t diff --git a/sac/src/vacusrpar.t.default b/sac/src/vacusrpar.t.default new file mode 100644 index 0000000..01275e2 --- /dev/null +++ b/sac/src/vacusrpar.t.default @@ -0,0 +1,3 @@ +! Include the gravity parameters + +INCLUDE:vacusrpar.gravity.t diff --git a/sac/vac.par b/sac/vac.par deleted file mode 120000 index 4a96520..0000000 --- a/sac/vac.par +++ /dev/null @@ -1 +0,0 @@ -par/mhdmodes \ No newline at end of file diff --git a/sac/vacini b/sac/vacini new file mode 120000 index 0000000..681bea7 --- /dev/null +++ b/sac/vacini @@ -0,0 +1 @@ +src/vacini \ No newline at end of file