diff --git a/examples/advection_rusanov/GNUmakefile b/examples/advection_rusanov/GNUmakefile new file mode 100644 index 00000000..691177ff --- /dev/null +++ b/examples/advection_rusanov/GNUmakefile @@ -0,0 +1,50 @@ +# NGA location if not yet defined +NGA_HOME ?= ../.. + +# Compilation parameters +PRECISION = DOUBLE +USE_MPI = TRUE +USE_FFTW = FALSE +USE_HYPRE = FALSE +USE_LAPACK = TRUE +PROFILE = FALSE +DEBUG = TRUE +COMP = gnu +EXEBASE = nga + +# Directories that contain user-defined code +Udirs := src + +# Include user-defined sources +Upack += $(foreach dir, $(Udirs), $(wildcard $(dir)/Make.package)) +Ulocs += $(foreach dir, $(Udirs), $(wildcard $(dir))) +include $(Upack) +INCLUDE_LOCATIONS += $(Ulocs) +VPATH_LOCATIONS += $(Ulocs) + +# Define external libraries - this can also be in .profile +HDF5_DIR = /opt/hdf5 +HYPRE_DIR = /opt/hypre +LAPACK_DIR = /opt/lapack-3.10.1 +FFTW_DIR = /opt/fftw-3.3.10 + +# NGA compilation definitions +include $(NGA_HOME)/tools/GNUMake/Make.defs + +# Include NGA base code +Bdirs := particles core hyperbolic data solver config grid libraries +Bpack += $(foreach dir, $(Bdirs), $(NGA_HOME)/src/$(dir)/Make.package) +include $(Bpack) + +# Inform user of Make.packages used +ifdef Ulocs + $(info Taking user code from: $(Ulocs)) +endif +$(info Taking base code from: $(Bdirs)) + +# Target definition +all: $(executable) + @echo COMPILATION SUCCESSFUL + +# NGA compilation rules +include $(NGA_HOME)/tools/GNUMake/Make.rules diff --git a/examples/advection_rusanov/README b/examples/advection_rusanov/README new file mode 100644 index 00000000..fb741241 --- /dev/null +++ b/examples/advection_rusanov/README @@ -0,0 +1,2 @@ +standard advection test case advecting different shapes over a periodic domain +to evaluate accuracy and diffusion diff --git a/examples/advection_rusanov/input b/examples/advection_rusanov/input new file mode 100644 index 00000000..267b7e03 --- /dev/null +++ b/examples/advection_rusanov/input @@ -0,0 +1,24 @@ +# Parallelization +Partition : 2 1 1 + +# Mesh definition +Lx : 1.0 +nx : 128 +Ly : 1.0 +ny : 128 +Lz : 1.0 +nz : 1 +Periodic : true + +# Initialization +Shapes : c + +# Direction +Advection direction : 1 1 0 + +# Time integration +Max cfl number : 0.5 + +# Ensight output +Ensight output period : 0.01 + diff --git a/examples/advection_rusanov/src/Make.package b/examples/advection_rusanov/src/Make.package new file mode 100644 index 00000000..a7a92785 --- /dev/null +++ b/examples/advection_rusanov/src/Make.package @@ -0,0 +1,2 @@ +# List here the extra files here +f90EXE_sources += geometry.f90 simulation.f90 diff --git a/examples/advection_rusanov/src/geometry.f90 b/examples/advection_rusanov/src/geometry.f90 new file mode 100644 index 00000000..b38098e3 --- /dev/null +++ b/examples/advection_rusanov/src/geometry.f90 @@ -0,0 +1,80 @@ +!> Various definitions and tools for initializing NGA2 config +module geometry + use config_class, only: config + use precision, only: WP + implicit none + private + + !> Single config + type(config), public :: cfg + + public :: geometry_init + +contains + + !> Initialization of problem geometry + subroutine geometry_init + use sgrid_class, only: sgrid + use param, only: param_read + implicit none + type(sgrid) :: grid + + ! Create a grid from input params + create_grid: block + use sgrid_class, only: cartesian + logical :: periodic + integer :: i, j, k, nx, ny, nz + real(WP) :: Lx, Ly, Lz + real(WP), dimension(:), allocatable :: x, y, z + + ! read in grid definition + call param_read('Lx', Lx) + call param_read('nx', nx) + call param_read('Ly', Ly) + call param_read('ny', ny) + call param_read('Lz', Lz) + call param_read('nz', nz) + call param_read('Periodic', periodic) + + ! allocate + allocate(x(nx+1), y(ny+1), z(nz+1)) + + ! create simple rectilinear grid + do i = 1, nx+1 + x(i) = real(i-1,WP) / real(nx, WP) * Lx + end do + do j = 1, ny+1 + y(j) = real(j-1,WP) / real(ny, WP) * Ly + end do + do k = 1, nz+1 + z(k) = real(k-1,WP) / real(nz, WP) * Lz + end do + + ! generate serial grid object + grid=sgrid(coord=cartesian, no=1, x=x, y=x, z=x, xper=periodic, & + & yper=periodic, zper=periodic) + + end block create_grid + + ! create a config from that grid on our entire group + create_cfg: block + use parallel, only: group + integer, dimension(3) :: partition + + ! read in partition + call param_read('Partition', partition, short='p') + + ! create partitioned grid + cfg = config(grp=group, decomp=partition, grid=grid) + + end block create_cfg + + ! Create masks for this config + create_walls: block + cfg%VF=1.0_WP + end block create_walls + + end subroutine geometry_init + +end module geometry + diff --git a/examples/advection_rusanov/src/simulation.f90 b/examples/advection_rusanov/src/simulation.f90 new file mode 100644 index 00000000..91435906 --- /dev/null +++ b/examples/advection_rusanov/src/simulation.f90 @@ -0,0 +1,310 @@ +!> Various definitions and tools for running an NGA2 simulation +module simulation + use precision, only: WP + use geometry, only: cfg + use rusanov_class, only: rusanov + use hyperbolic_advection, only: make_advec_rusanov + use timetracker_class, only: timetracker + use ensight_class, only: ensight + use partmesh_class, only: partmesh + use event_class, only: event + use monitor_class, only: monitor + implicit none + private + + !> Single-phase incompressible flow solver, pressure and implicit solvers, and a time tracker + type(rusanov), public :: fs + type(timetracker), public :: time + + !> Ensight postprocessing + type(ensight) :: ens_out + type(event) :: ens_evt + + !> simulation monitor file + type(monitor) :: mfile, cflfile + + !> simulation control functions + public :: simulation_init, simulation_run, simulation_final + +contains + + !> circles and donuts shape + ! TODO fill between zero and one on boundaries + subroutine init_donut(iR, oR, U) + implicit none + real(WP), dimension(cfg%imino_:,cfg%jmino_:,cfg%kmino_:), intent(out) :: U + real(WP), intent(in) :: iR, oR + real(WP), dimension(3) :: c, x + real(WP) :: r2 + integer :: i, j, k + + U(:,:,:) = 0.0_WP + + c = 0.5 * (/ cfg%xL, cfg%yL, cfg%zL /) + + do k = cfg%kmin_, cfg%kmax_ + x(3) = cfg%zm(k) + do j = cfg%jmin_, cfg%jmax_ + x(2) = cfg%ym(j) + do i = cfg%imin_, cfg%imax_ + x(1) = cfg%xm(i) + r2 = sum((c - x)**2) + if (iR**2 .lt. r2 .and. r2 .lt. oR**2) U(i,j,k) = 1.0_WP + end do + end do + end do + + end subroutine init_donut + + !> bar so we can test 1d behavior + subroutine init_bar(xmin, xmax, U) + implicit none + real(WP), dimension(cfg%imino_:,cfg%jmino_:,cfg%kmino_:), intent(out) :: U + real(WP), intent(in) :: xmin, xmax + integer :: i + + U(:,:,:) = 0.0_WP + + do i = cfg%imin_, cfg%imax_ + if (xmin .lt. cfg%xm(i) .and. cfg%xm(i) .lt. xmax) then + U(i,:,:) = 1.0_WP + end if + end do + + end subroutine + + !> initialization of problem solver + subroutine simulation_init() + use param, only: param_read + use string, only: str_medium + use messager, only: die + implicit none + + integer :: numfields + character, dimension(:), allocatable :: fields + + ! read shapes + read_shapes: block + character(len=str_medium) :: fieldbuffer + integer :: i, j + + call param_read('Shapes', fieldbuffer) + fieldbuffer = adjustl(fieldbuffer) + + ! check for duplicates + do i = 1, str_medium + if (fieldbuffer(i:i) .ne. ' ') then + do j = i+1, str_medium + if (fieldbuffer(i:i) .eq. fieldbuffer(j:j)) then + call die("field name '" // fieldbuffer(i:i) // "' appears twice") + end if + end do + end if + end do + + ! count valid fields and warn user + numfields = 0 + do i = 1, str_medium + select case (fieldbuffer(i:i)) + case ('c') ! circle + numfields = numfields + 1 + case ('d') ! donut + numfields = numfields + 1 + case ('b') ! bar + numfields = numfields + 1 + case (' ') ! skip spaces + case default + write(*, *) "unrecognized shape: '", fieldbuffer(i:i), "'" + end select + end do + if (numfields .eq. 0) call die("must advect at least one shape") + + ! read fields into array + allocate(fields(numfields)) + numfields = 0 + do i = 1, str_medium + select case (fieldbuffer(i:i)) + case ('c') ! circle + numfields = numfields + 1 + fields(numfields) = 'c' + case ('d') ! donut + numfields = numfields + 1 + fields(numfields) = 'd' + case ('b') ! bar + numfields = numfields + 1 + fields(numfields) = 'b' + case default + end select + end do + + end block read_shapes + + !TODO need to figure out what the right interaction with this is + initialize_timetracker: block + + time = timetracker(amRoot=cfg%amRoot) + call param_read('Max cfl number', time%cflmax) + time%dt = 1e-3_WP + time%itmax = 1 + + end block initialize_timetracker + + ! create a single-phase flow solver without bconds + create_and_initialize_flow_solver: block + real(WP), dimension(3) :: vel + real(WP) :: traversal_length + + ! read velocity direction, find velocity + call param_read('Advection direction', vel, 'Advec dir') + traversal_length = sqrt(sum(vel**2)) * cfg%xL + vel = vel * traversal_length + + ! call constructor + fs = make_advec_rusanov(cfg, numfields, vel) + + end block create_and_initialize_flow_solver + + ! prepare initial fields + initialize_fields: block + integer :: i + real(WP), dimension(:,:,:), pointer :: fieldptr + + do i = 1, numfields + fieldptr => fs%Uc(i,:,:,:) + select case (fields(i)) + case ('c') ! circle + call init_donut(0.0_WP, 0.375_WP * cfg%xL, fieldptr) + case ('d') ! donut + call init_donut(0.250_WP * cfg%xL, 0.375_WP * cfg%xL, fieldptr) + case ('b') ! bar + call init_bar(0.333_WP * cfg%xL, 0.666_WP * cfg%xL, fieldptr) + end select + end do + + end block initialize_fields + + ! Add Ensight output + create_ensight: block + use string, only: str_short + integer :: i + character(len=4) :: istr + real(WP), dimension(:,:,:), pointer :: ptr1, ptr2, ptr3 + + ! Create Ensight output from cfg + ens_out=ensight(cfg=cfg, name='AdvectionTest') + + ! Create event for Ensight output + ens_evt=event(time=time, name='Ensight output') + call param_read('Ensight output period', ens_evt%tper) + + ! Add variables to output + do i = 1, numfields + ptr1 => fs%Uc(i,:,:,:) + write(istr,'(I2.2)') i + call ens_out%add_scalar('U'//istr, ptr1) + end do + + ptr1 => fs%params(1,:,:,:) + ptr2 => fs%params(2,:,:,:) + ptr3 => fs%params(3,:,:,:) + call ens_out%add_vector('velocity', ptr1, ptr2, ptr3) + + ! Output to ensight + if (ens_evt%occurs()) call ens_out%write_data(time%t) + + end block create_ensight + + ! Create a monitor file + create_monitor: block + use string, only: str_short + character(len=str_short) :: fieldname + integer :: i + real(WP), pointer :: real_ptr + + ! Prepare some info about fields + call fs%get_cfl(time%dt,time%cfl) + call fs%get_range() + + ! Create simulation monitor + mfile = monitor(fs%cfg%amRoot, 'simulation') + call mfile%add_column(time%n, 'Timestep number') + call mfile%add_column(time%t, 'Time') + call mfile%add_column(time%dt, 'Timestep size') + call mfile%add_column(time%cfl, 'Maximum CFL') + do i = 1, numfields + write(fieldname,'("[",a,"]min")') fields(i) + real_ptr => fs%Umin(i) + call mfile%add_column(real_ptr, fieldname) + write(fieldname,'("[",a,"]max")') fields(i) + real_ptr => fs%Umax(i) + call mfile%add_column(real_ptr, fieldname) + end do + call mfile%write() + + ! Create CFL monitor + cflfile = monitor(fs%cfg%amRoot, 'cfl') + call cflfile%add_column(time%n, 'Timestep number') + call cflfile%add_column(time%t, 'Time') + call cflfile%add_column(fs%CFL_x, 'CFLx') + call cflfile%add_column(fs%CFL_y, 'CFLy') + call cflfile%add_column(fs%CFL_z, 'CFLz') + call cflfile%write() + + end block create_monitor + + end subroutine simulation_init + + !> do the thing + subroutine simulation_run + implicit none + + ! Perform time integration + do while (.not.time%done()) + + ! Increment time + call fs%get_cfl(time%dt, time%cfl) + call time%adjust_dt() + call time%increment() + + ! take step (Strang) + fs%dU(:, :, :, :) = 0.0_WP + !call fs%apply_bcond(time%t, time%dt) + call fs%calc_dU_x(0.5 * time%dt) + fs%Uc = fs%Uc + fs%dU + fs%dU(:, :, :, :) = 0.0_WP + !call fs%apply_bcond(time%t, time%dt) + call fs%calc_dU_y(time%dt) + fs%Uc = fs%Uc + fs%dU + fs%dU(:, :, :, :) = 0.0_WP + !call fs%apply_bcond(time%t, time%dt) + call fs%calc_dU_x(0.5 * time%dt) + fs%Uc = fs%Uc + fs%dU + + ! Output to ensight + if (ens_evt%occurs()) call ens_out%write_data(time%t) + + ! Perform and output monitoring + call fs%get_range() + call mfile%write() + call cflfile%write() + + end do + + end subroutine simulation_run + + !> Finalize the NGA2 simulation + subroutine simulation_final + implicit none + + ! Get rid of all objects - need destructors + !TODO + ! monitor + ! ensight + ! bcond + ! timetracker + + + end subroutine simulation_final + +end module simulation + diff --git a/examples/sod_rusanov/GNUmakefile b/examples/sod_rusanov/GNUmakefile new file mode 100644 index 00000000..691177ff --- /dev/null +++ b/examples/sod_rusanov/GNUmakefile @@ -0,0 +1,50 @@ +# NGA location if not yet defined +NGA_HOME ?= ../.. + +# Compilation parameters +PRECISION = DOUBLE +USE_MPI = TRUE +USE_FFTW = FALSE +USE_HYPRE = FALSE +USE_LAPACK = TRUE +PROFILE = FALSE +DEBUG = TRUE +COMP = gnu +EXEBASE = nga + +# Directories that contain user-defined code +Udirs := src + +# Include user-defined sources +Upack += $(foreach dir, $(Udirs), $(wildcard $(dir)/Make.package)) +Ulocs += $(foreach dir, $(Udirs), $(wildcard $(dir))) +include $(Upack) +INCLUDE_LOCATIONS += $(Ulocs) +VPATH_LOCATIONS += $(Ulocs) + +# Define external libraries - this can also be in .profile +HDF5_DIR = /opt/hdf5 +HYPRE_DIR = /opt/hypre +LAPACK_DIR = /opt/lapack-3.10.1 +FFTW_DIR = /opt/fftw-3.3.10 + +# NGA compilation definitions +include $(NGA_HOME)/tools/GNUMake/Make.defs + +# Include NGA base code +Bdirs := particles core hyperbolic data solver config grid libraries +Bpack += $(foreach dir, $(Bdirs), $(NGA_HOME)/src/$(dir)/Make.package) +include $(Bpack) + +# Inform user of Make.packages used +ifdef Ulocs + $(info Taking user code from: $(Ulocs)) +endif +$(info Taking base code from: $(Bdirs)) + +# Target definition +all: $(executable) + @echo COMPILATION SUCCESSFUL + +# NGA compilation rules +include $(NGA_HOME)/tools/GNUMake/Make.rules diff --git a/examples/sod_rusanov/input_aligned b/examples/sod_rusanov/input_aligned new file mode 100644 index 00000000..4ff2de33 --- /dev/null +++ b/examples/sod_rusanov/input_aligned @@ -0,0 +1,21 @@ +# Parallelization +Partition : 2 1 1 + +# Mesh definition +Lx : 1.0 +nx : 64 +Ly : 1.0 +ny : 64 + +# Mesh Angle +Mesh angle : 0.0 + +# Initialization +Initial x velocity : 0.0 + +# Time integration +Max cfl number : 0.98 + +# Ensight output +Ensight output period : 0.01 + diff --git a/examples/sod_rusanov/input_skew b/examples/sod_rusanov/input_skew new file mode 100644 index 00000000..5974f5f3 --- /dev/null +++ b/examples/sod_rusanov/input_skew @@ -0,0 +1,21 @@ +# Parallelization +Partition : 2 2 1 + +# Mesh definition +Lx : 1.0 +nx : 64 +Ly : 1.0 +ny : 64 + +# Mesh Angle +Mesh angle : 35.0 + +# Initialization +Initial x velocity : 0.0 + +# Time integration +Max cfl number : 0.98 + +# Ensight output +Ensight output period : 0.01 + diff --git a/examples/sod_rusanov/src/Make.package b/examples/sod_rusanov/src/Make.package new file mode 100644 index 00000000..a7a92785 --- /dev/null +++ b/examples/sod_rusanov/src/Make.package @@ -0,0 +1,2 @@ +# List here the extra files here +f90EXE_sources += geometry.f90 simulation.f90 diff --git a/examples/sod_rusanov/src/geometry.f90 b/examples/sod_rusanov/src/geometry.f90 new file mode 100644 index 00000000..6a93a5c5 --- /dev/null +++ b/examples/sod_rusanov/src/geometry.f90 @@ -0,0 +1,80 @@ +!> Various definitions and tools for initializing NGA2 config +module geometry + use config_class, only: config + use precision, only: WP + implicit none + private + + !> Single config + type(config), public :: cfg + + public :: geometry_init + +contains + + !> Initialization of problem geometry + subroutine geometry_init + use sgrid_class, only: sgrid + use param, only: param_read + implicit none + type(sgrid) :: grid + + ! Create a grid from input params + create_grid: block + use sgrid_class, only: cartesian + integer :: i, j, k, nx, ny, nz + real(WP) :: Lx, Ly, Lz + real(WP), dimension(:), allocatable :: x, y, z + + ! read in grid definition + call param_read('Lx', Lx) + call param_read('nx', nx) + call param_read('Ly', Ly) + call param_read('ny', ny) + nz = 1 + Lz = 1.0_WP + + ! allocate + allocate(x(nx+1), y(ny+1), z(nz+1)) + + ! create simple rectilinear grid + do i = 1, nx+1 + x(i) = real(i-1-nx/2,WP) / real(nx, WP) * Lx + end do + do j = 1, ny+1 + y(j) = real(j-1-nz/2,WP) / real(ny, WP) * Ly + end do + do k = 1, nz+1 + z(k) = real(k-1,WP) / real(nz, WP) * Lz + end do + + ! generate serial grid object + grid=sgrid(coord=cartesian, no=1, x=x, y=x, z=x, xper=.false., & + & yper=.false., zper=.false.) + + deallocate(x, y, z) + + end block create_grid + + ! create a config from that grid on our entire group + create_cfg: block + use parallel, only: group + integer, dimension(3) :: partition + + ! read in partition + call param_read('Partition', partition, short='p') + + ! create partitioned grid + cfg = config(grp=group, decomp=partition, grid=grid) + + end block create_cfg + + ! Create masks for this config + create_walls: block + cfg%VF = 1.0_WP + end block create_walls + + end subroutine geometry_init + +end module geometry + diff --git a/examples/sod_rusanov/src/simulation.f90 b/examples/sod_rusanov/src/simulation.f90 new file mode 100644 index 00000000..c0d58da3 --- /dev/null +++ b/examples/sod_rusanov/src/simulation.f90 @@ -0,0 +1,333 @@ +!> Various definitions and tools for running an NGA2 simulation +module simulation + use precision, only: WP + use geometry, only: cfg + use rusanov_class, only: rusanov, EXTRAPOLATE + use hyperbolic_euler, only: make_euler_rusanov, euler_tocons, & + & euler_tophys, SOD_PHYS_L, SOD_PHYS_R, DIATOMIC_GAMMA + use timetracker_class, only: timetracker + use ensight_class, only: ensight + use partmesh_class, only: partmesh + use event_class, only: event + use monitor_class, only: monitor + use string, only: str_medium + implicit none + private + + !> Flow solver and a time tracker + type(rusanov), public :: fs + type(timetracker), public :: time + + !> Ensight postprocessing + character(len=str_medium) :: ens_out_name + type(ensight) :: ens_out + type(event) :: ens_evt + + !> physical value arrays for ensight output + real(WP), dimension(:,:,:,:), pointer :: phys_out + + !> simulation monitor file + type(monitor) :: mfile, cflfile, consfile + + !> simulation control functions + public :: simulation_init, simulation_run, simulation_final + +contains + + !> update phys + subroutine update_phys_out() + implicit none + integer :: i, j, k + + do k = cfg%kmin_, cfg%kmax_ + do j = cfg%jmin_, cfg%jmax_ + do i = cfg%imin_, cfg%imax_ + call euler_tophys(DIATOMIC_GAMMA, fs%Uc(:,i,j,k), phys_out(1:5,i,j,k)) + phys_out(6,i,j,k) = sqrt(sum(phys_out(2:4,i,j,k)**2) / (DIATOMIC_GAMMA * phys_out(5,i,j,k) / phys_out(1,i,j,k))) + end do + end do + end do + + call cfg%sync(phys_out) + + end subroutine + + !> functions localize the sides of the domain + function side_locator_x_l(pg, i, j, k) result(loc) + use pgrid_class, only: pgrid + implicit none + class(pgrid), intent(in) :: pg + integer, intent(in) :: i, j, k + logical :: loc + loc = .false. + if (i.lt.pg%imin) loc = .true. + end function side_locator_x_l + function side_locator_x_r(pg, i, j, k) result(loc) + use pgrid_class, only: pgrid + implicit none + class(pgrid), intent(in) :: pg + integer, intent(in) :: i, j, k + logical :: loc + loc = .false. + if (i.gt.pg%imax) loc = .true. + end function side_locator_x_r + function side_locator_y_l(pg, i, j, k) result(loc) + use pgrid_class, only: pgrid + implicit none + class(pgrid), intent(in) :: pg + integer, intent(in) :: i, j, k + logical :: loc + loc = .false. + if (j.lt.pg%jmin) loc = .true. + end function side_locator_y_l + function side_locator_y_r(pg, i, j, k) result(loc) + use pgrid_class, only: pgrid + implicit none + class(pgrid), intent(in) :: pg + integer, intent(in) :: i, j, k + logical :: loc + loc = .false. + if (j.gt.pg%jmax) loc = .true. + end function side_locator_y_r + + !> initialization of problem solver + subroutine simulation_init() + use param, only: param_read + use messager, only: die + implicit none + + !TODO need to figure out what the right interaction with this is + initialize_timetracker: block + + time = timetracker(amRoot=cfg%amRoot) + call param_read('Max cfl number', time%cflmax) + time%dt = 1e-3_WP + time%itmax = 1 + + end block initialize_timetracker + + ! create a single-phase flow solver + create_and_initialize_flow_solver: block + + ! call constructor + fs = make_euler_rusanov(cfg) + + ! set bcs + call fs%add_bcond(name='openxl' , type=EXTRAPOLATE, & + & locator=side_locator_x_l, dir='xl') + call fs%add_bcond(name='openxr' , type=EXTRAPOLATE, & + & locator=side_locator_x_r, dir='xr') + call fs%add_bcond(name='openyl' , type=EXTRAPOLATE, & + & locator=side_locator_y_l, dir='yl') + call fs%add_bcond(name='openyr' , type=EXTRAPOLATE, & + & locator=side_locator_y_r, dir='yr') + + end block create_and_initialize_flow_solver + + ! prepare initial fields + initialize_fields: block + real(WP), dimension(5) :: lval, rval, val + real(WP) :: angle, x0, y0, x1, y1, w + integer :: i, j, k + + call param_read('Mesh angle', angle) + + write(ens_out_name,'(A,F0.1)') 'SodAngle', angle + + angle = angle * 2.0_WP * 4.0_WP * atan(1.0_WP) / 360 + + call euler_tocons(DIATOMIC_GAMMA, SOD_PHYS_L, lval) + call euler_tocons(DIATOMIC_GAMMA, SOD_PHYS_R, rval) + + do j = cfg%jmino_, cfg%jmaxo_ + y0 = cfg%y(j); y1 = cfg%y(j+1); + do i = cfg%imino_, cfg%imaxo_ + x0 = cfg%x(i); x1 = cfg%x(i+1); + if (cos(angle) * x1 + sin(angle) * y1 .le. 0.0_WP) then + val = lval + else if (cos(angle) * x0 + sin(angle) * y0 .ge. 0.0_WP) then + val = rval + else + w = 0.5_WP * (cotan(angle) * x0 + y0) * (tan(angle) * y0 + x0) + if (-cotan(angle) * x0 .gt. y1) then + w = w - 0.5_WP * (tan(angle) * y1 + x0) * (cotan(angle) * x0 + y1) + end if + if (-tan(angle) * y0 .gt. x1) then + w = w - 0.5_WP * (tan(angle) * y0 + x1) * (cotan(angle) * x1 + y0) + end if + w = w / (cfg%dx(i) * cfg%dy(j)) + val = lval * w + rval * (1.0_WP - w) + end if + do k = cfg%kmino_, cfg%kmaxo_ + fs%Uc(:,i,j,k) = val + end do + end do + end do + + end block initialize_fields + + ! Add Ensight output + create_ensight: block + use string, only: str_short + real(WP), dimension(:,:,:), pointer :: ptr1, ptr2, ptr3 + + ! create array to hold physical coordinates + allocate(phys_out(6,cfg%imino_:cfg%imaxo_,cfg%jmino_:cfg%jmaxo_, & + & cfg%kmino_:cfg%kmaxo_)) + + ! Create Ensight output from cfg + ens_out = ensight(cfg=cfg, name=ens_out_name) + + ! Create event for Ensight output + ens_evt = event(time=time, name='Ensight output') + call param_read('Ensight output period', ens_evt%tper) + + ! Add variables to output + ptr1 => phys_out(1,:,:,:) + call ens_out%add_scalar('density', ptr1) + ptr1 => phys_out(2,:,:,:) + ptr2 => phys_out(3,:,:,:) + ptr3 => phys_out(4,:,:,:) + call ens_out%add_vector('velocity', ptr1, ptr2, ptr3) + ptr1 => phys_out(5,:,:,:) + call ens_out%add_scalar('pressure', ptr1) + ptr1 => fs%params(1,:,:,:) + call ens_out%add_scalar('gamma', ptr1) + ptr1 => phys_out(6,:,:,:) + call ens_out%add_scalar('Ma', ptr1) + + ! Output to ensight + if (ens_evt%occurs()) then + call update_phys_out() + call ens_out%write_data(time%t) + end if + + end block create_ensight + + ! Create a monitor file + create_monitor: block + use string, only: str_short + real(WP), pointer :: real_ptr + + ! Prepare some info about fields + call fs%get_cfl(time%dt, time%cfl) + call fs%get_range() + + ! Create simulation monitor + mfile = monitor(fs%cfg%amRoot, 'simulation') + call mfile%add_column(time%n, 'Timestep number') + call mfile%add_column(time%t, 'Time') + call mfile%add_column(time%dt, 'Timestep size') + call mfile%add_column(time%cfl, 'Maximum CFL') + ! not sure why this doesn't work? + !call mfile%add_column(real_ptr, fields(i:i)//'min') + !call mfile%add_column(real_ptr, fields(i:i)//'max') + real_ptr => fs%Umin(1) + call mfile%add_column(real_ptr, 'dens_min') + real_ptr => fs%Umax(1) + call mfile%add_column(real_ptr, 'dens_max') + real_ptr => fs%Umin(2) + call mfile%add_column(real_ptr, 'momx_min') + real_ptr => fs%Umax(2) + call mfile%add_column(real_ptr, 'momx_max') + real_ptr => fs%Umin(3) + call mfile%add_column(real_ptr, 'momy_min') + real_ptr => fs%Umax(3) + call mfile%add_column(real_ptr, 'momy_max') + real_ptr => fs%Umin(4) + call mfile%add_column(real_ptr, 'momz_min') + real_ptr => fs%Umax(4) + call mfile%add_column(real_ptr, 'momz_max') + real_ptr => fs%Umin(5) + call mfile%add_column(real_ptr, 'totE_min') + real_ptr => fs%Umax(5) + call mfile%add_column(real_ptr, 'totE_max') + call mfile%write() + + ! Create CFL monitor + cflfile = monitor(fs%cfg%amRoot, 'cfl') + call cflfile%add_column(time%n, 'Timestep number') + call cflfile%add_column(time%t, 'Time') + call cflfile%add_column(fs%CFL_x, 'CFLx') + call cflfile%add_column(fs%CFL_y, 'CFLy') + call cflfile%add_column(fs%CFL_z, 'CFLz') + call cflfile%write() + + ! Create conservation monitor + consfile = monitor(fs%cfg%amRoot, 'conservation') + call consfile%add_column(time%n, 'Timestep number') + call consfile%add_column(time%t, 'Time') + real_ptr => fs%Uint(1) + call consfile%add_column(real_ptr, 'dens_int') + real_ptr => fs%Uint(2) + call consfile%add_column(real_ptr, 'momx_int') + real_ptr => fs%Uint(3) + call consfile%add_column(real_ptr, 'momy_int') + real_ptr => fs%Uint(4) + call consfile%add_column(real_ptr, 'momz_int') + real_ptr => fs%Uint(5) + call consfile%add_column(real_ptr, 'totE_int') + call consfile%write() + + end block create_monitor + + end subroutine simulation_init + + !> do the thing + subroutine simulation_run + implicit none + + ! Perform time integration + do while (.not.time%done()) + + ! Increment time + call fs%get_cfl(time%dt, time%cfl) + call time%adjust_dt() + call time%increment() + + ! take step (Strang) + call fs%apply_bcond(time%t, time%dt) + fs%dU(:, :, :, :) = 0.0_WP + call fs%calc_dU_x(0.5 * time%dt) + fs%Uc = fs%Uc + fs%dU + call fs%apply_bcond(time%t, time%dt) + fs%dU(:, :, :, :) = 0.0_WP + call fs%calc_dU_y(time%dt) + fs%Uc = fs%Uc + fs%dU + call fs%apply_bcond(time%t, time%dt) + fs%dU(:, :, :, :) = 0.0_WP + call fs%calc_dU_x(0.5 * time%dt) + fs%Uc = fs%Uc + fs%dU + + ! Output to ensight + if (ens_evt%occurs()) then + call update_phys_out() + call ens_out%write_data(time%t) + end if + + ! Perform and output monitoring + call fs%get_range() + call mfile%write() + call cflfile%write() + call consfile%write() + + end do + + end subroutine simulation_run + + !> Finalize the NGA2 simulation + subroutine simulation_final + implicit none + + ! Get rid of all objects - need destructors + !TODO + ! monitor + ! ensight + ! bcond + ! timetracker + + + end subroutine simulation_final + +end module simulation + diff --git a/src/hyperbolic/Make.package b/src/hyperbolic/Make.package index 73e3d8da..379b832a 100644 --- a/src/hyperbolic/Make.package +++ b/src/hyperbolic/Make.package @@ -1,5 +1,5 @@ -f90EXE_sources += hyperbolic_general.f90 muscl_class.f90 advection.f90 euler.f90 -# houssem2d.f90 rusanov_class.f90 +f90EXE_sources += hyperbolic_general.f90 muscl_class.f90 advection.f90 euler.f90 rusanov_class.f90 +# houssem2d.f90 INCLUDE_LOCATIONS += $(NGA_HOME)/src/hyperbolic VPATH_LOCATIONS += $(NGA_HOME)/src/hyperbolic diff --git a/src/hyperbolic/advection.f90 b/src/hyperbolic/advection.f90 index 0848f1a9..59353540 100644 --- a/src/hyperbolic/advection.f90 +++ b/src/hyperbolic/advection.f90 @@ -5,20 +5,22 @@ !> Originally written by John P Wakefield in December 2022. module hyperbolic_advection - use precision, only: WP - use string, only: str_medium - use config_class, only: config - use hyperbolic, only: eigenvals_ftype, rsolver_ftype, limiter_ftype - use muscl_class, only: muscl + use precision, only: WP + use string, only: str_medium + use config_class, only: config + use hyperbolic, only: eigenvals_ftype, rsolver_ftype, limiter_ftype, flux_ftype + use muscl_class, only: muscl + use rusanov_class, only: rusanov implicit none real(WP), parameter :: advec_muscl_cflsafety = 0.99_WP real(WP), parameter :: advec_muscl_divzero_eps = 1e-12_WP - character(len=str_medium), parameter :: advec_muscl_name = 'MUSCL_CONST_ADVEC' + character(len=str_medium), parameter :: ADVEC_MUSCL_NAME = 'MUSCL_CONST_ADVEC' + character(len=str_medium), parameter :: ADVEC_RUS_NAME = 'RUSANOV_CONST_ADVEC' contains - !> factory + !> muscl factory function make_advec_muscl(cfg, N, limiter, velocity) result(solver) implicit none type(muscl) :: solver @@ -31,7 +33,7 @@ function make_advec_muscl(cfg, N, limiter, velocity) result(solver) procedure(rsolver_ftype), pointer :: rsolv_x_ptr, rsolv_y_ptr, rsolv_z_ptr integer :: i - name_actual = advec_muscl_name + name_actual = ADVEC_MUSCL_NAME evals_x_ptr => advec_evals_x; rsolv_x_ptr => advec_rsolv_x; evals_y_ptr => advec_evals_y; rsolv_y_ptr => advec_rsolv_y; @@ -54,6 +56,37 @@ function make_advec_muscl(cfg, N, limiter, velocity) result(solver) end function make_advec_muscl + !> rusanov factory + function make_advec_rusanov(cfg, N, velocity) result(solver) + implicit none + type(rusanov) :: solver + class(config), target, intent(in) :: cfg + integer, intent(in) :: N + real(WP), dimension(3), intent(in) :: velocity + procedure(eigenvals_ftype), pointer :: evals_x_ptr, evals_y_ptr, evals_z_ptr + procedure(flux_ftype), pointer :: flux_x_ptr, flux_y_ptr, flux_z_ptr + integer :: i + + evals_x_ptr => advec_evals_x; flux_x_ptr => advec_flux_x; + evals_y_ptr => advec_evals_y; flux_y_ptr => advec_flux_y; + evals_z_ptr => advec_evals_z; flux_z_ptr => advec_flux_z; + + ! build solver + solver = rusanov(cfg, ADVEC_RUS_NAME, N, 3, evals_x_ptr, evals_y_ptr, & + & evals_z_ptr, flux_x_ptr, flux_y_ptr, flux_z_ptr) + + ! put velocity in param array + do i = 1, 3 + solver%params(i,:,:,:) = velocity(i) + end do + + ! set velocity mask + solver%vel_mask_x(:) = .false. + solver%vel_mask_y(:) = .false. + solver%vel_mask_z(:) = .false. + + end function make_advec_rusanov + pure subroutine advec_evals_x(P, N, params, U, evals) implicit none integer, intent(in) :: P, N @@ -96,6 +129,39 @@ pure subroutine advec_evals_z(P, N, params, U, evals) end subroutine advec_evals_z + pure subroutine advec_flux_x(P, N, params, U, flux) + implicit none + integer, intent(in) :: P, N + real(WP), dimension(P), intent(in) :: params + real(WP), dimension(N), intent(in) :: u + real(WP), dimension(N), intent(out) :: flux + + flux(:) = params(1) * U(:) + + end subroutine advec_flux_x + + pure subroutine advec_flux_y(P, N, params, U, flux) + implicit none + integer, intent(in) :: P, N + real(WP), dimension(P), intent(in) :: params + real(WP), dimension(N), intent(in) :: u + real(WP), dimension(N), intent(out) :: flux + + flux(:) = params(2) * U(:) + + end subroutine advec_flux_y + + pure subroutine advec_flux_z(P, N, params, U, flux) + implicit none + integer, intent(in) :: P, N + real(WP), dimension(P), intent(in) :: params + real(WP), dimension(N), intent(in) :: u + real(WP), dimension(N), intent(out) :: flux + + flux(:) = params(3) * U(:) + + end subroutine advec_flux_z + pure subroutine advec_rsolv_simple(v, Ul, Ur, rs) real(WP), intent(in) :: v real(WP), dimension(:), intent(in) :: Ul, Ur diff --git a/src/hyperbolic/euler.f90 b/src/hyperbolic/euler.f90 index b281eb92..7b76caa5 100644 --- a/src/hyperbolic/euler.f90 +++ b/src/hyperbolic/euler.f90 @@ -10,7 +10,7 @@ module hyperbolic_euler use config_class, only: config use hyperbolic, only: VANLEER, eigenvals_ftype, rsolver_ftype, limiter_ftype, flux_ftype use muscl_class, only: muscl - !use rusanov_class, only: rusanov + use rusanov_class, only: rusanov implicit none real(WP), parameter :: euler_muscl_cflsafety = 0.92_WP @@ -72,41 +72,38 @@ function make_euler_muscl(cfg, limiter, gma) result(solver) end function make_euler_muscl !> rusanov factory - !function make_euler_rusanov(cfg, gma) result(solver) - ! implicit none - ! type(rusanov) :: solver - ! class(config), target, intent(in) :: cfg - ! real(WP), optional, intent(in) :: gma - ! real(WP) :: gma_actual - ! procedure(eigenvals_ftype), pointer :: evals_x_ptr, evals_y_ptr, evals_z_ptr - ! procedure(flux_ftype), pointer :: flux_x_ptr, flux_y_ptr, flux_z_ptr - - ! if (present(gma)) then - ! gma_actual = gma - ! else - ! gma_actual = DIATOMIC_GAMMA - ! end if - - ! evals_x_ptr => euler_evals_x - ! evals_y_ptr => euler_evals_y - ! evals_z_ptr => euler_evals_z - ! flux_x_ptr => euler_flux_x - ! flux_y_ptr => euler_flux_y - ! flux_z_ptr => euler_flux_z - - ! ! build solver - ! solver = rusanov(cfg, 'EULER_RUSANOV', 5, 1, evals_x_ptr, evals_y_ptr, & - ! & evals_z_ptr, flux_x_ptr, flux_y_ptr, flux_z_ptr) - - ! ! set param array to hold gamma - ! solver%params(1,:,:,:) = gma_actual - - ! ! set velocity mask - ! solver%vel_mask_x(:) = (/ .false., .true., .false., .false., .false. /) - ! solver%vel_mask_y(:) = (/ .false., .false., .true., .false., .false. /) - ! solver%vel_mask_z(:) = (/ .false., .false., .false., .true., .false. /) - - !end function make_euler_rusanov + function make_euler_rusanov(cfg, gma) result(solver) + implicit none + type(rusanov) :: solver + class(config), target, intent(in) :: cfg + real(WP), optional, intent(in) :: gma + real(WP) :: gma_actual + procedure(eigenvals_ftype), pointer :: evals_x_ptr, evals_y_ptr, evals_z_ptr + procedure(flux_ftype), pointer :: flux_x_ptr, flux_y_ptr, flux_z_ptr + + if (present(gma)) then + gma_actual = gma + else + gma_actual = DIATOMIC_GAMMA + end if + + evals_x_ptr => euler_evals_x; flux_x_ptr => euler_flux_x; + evals_y_ptr => euler_evals_y; flux_y_ptr => euler_flux_y; + evals_z_ptr => euler_evals_z; flux_z_ptr => euler_flux_z; + + ! build solver + solver = rusanov(cfg, 'EULER_RUSANOV', 5, 1, evals_x_ptr, evals_y_ptr, & + & evals_z_ptr, flux_x_ptr, flux_y_ptr, flux_z_ptr) + + ! set param array to hold gamma + solver%params(1,:,:,:) = gma_actual + + ! set velocity mask + solver%vel_mask_x(:) = (/ .false., .true., .false., .false., .false. /) + solver%vel_mask_y(:) = (/ .false., .false., .true., .false., .false. /) + solver%vel_mask_z(:) = (/ .false., .false., .false., .true., .false. /) + + end function make_euler_rusanov !> Convert to Physical Coordinates (velocity and pressure to momentum and energy) pure subroutine euler_tophys(gma, cons, phys) @@ -150,7 +147,8 @@ pure subroutine euler_flux_x(P, N, params, U, flux) real(WP), dimension(N), intent(out) :: flux real(WP) :: pressure - pressure = (params(1) - 1.0_WP) * (U(3) - 0.5_WP * U(2)**2 / U(1)) + pressure = (params(1) - 1.0_WP) * (U(5) - 0.5_WP * (U(2)**2 + U(3)**2 + & + U(4)**2)/ U(1)) flux(1) = U(2) flux(2) = U(2)**2 / U(1) + pressure diff --git a/src/hyperbolic/rusanov_class.f90 b/src/hyperbolic/rusanov_class.f90 new file mode 100644 index 00000000..d10857d9 --- /dev/null +++ b/src/hyperbolic/rusanov_class.f90 @@ -0,0 +1,561 @@ +!> Rusanov solver class +!> +!> This first-order Rusanov scheme is extremely diffusive, but foolproof (as +!> long as a flux function is known. +!> +!> Originally written by John P Wakefield, Feb 2023. +!> +!> +module rusanov_class + use precision, only: WP + use string, only: str_medium + use config_class, only: config + use iterator_class, only: iterator + use hyperbolic, only: eigenvals_ftype, flux_ftype + implicit none + private + + ! expose type/constructor/methods + public :: rusanov, rus_bc + + ! expose interfaces + public :: eigenvals_ftype, flux_ftype + + ! List of known available bcond types for this solver + integer(1), parameter, public :: EXTRAPOLATE = 1_1 !< 0th order extrapolation + integer(1), parameter, public :: REFLECT = 2_1 !< reflection by velmask + + !> boundary conditions for the hyperbolic solver + type :: rus_bc + type(rus_bc), pointer :: next !< linked list of bconds + character(len=str_medium) :: name = 'UNNAMED_BCOND' !< bcond name (default UNNAMED_BCOND) + integer(1) :: type !< bcond type + type(iterator) :: itr !< this is the iterator for the bcond + integer(1) :: dir !< bcond direction 0-6 + end type rus_bc + + !> solver object definition + type :: rusanov + + ! config / pgrid object + class(config), pointer :: cfg !< config / pgrid information + + ! name + character(len=str_medium) :: name = 'UNNAMED_RUSANOV' !< solver name + + ! system information + integer :: N !< system dimension + integer :: P !< number of parameters + procedure(eigenvals_ftype), pointer, nopass :: evals_x, evals_y, evals_z + procedure(flux_ftype), pointer, nopass :: flux_x, flux_y, flux_z + logical, dimension(:), allocatable :: vel_mask_x, vel_mask_y, vel_mask_z + + ! boundary condition list + integer :: nbc !< number of bcond for our solver + type(rus_bc), pointer :: first_bc !< list of bcond for our solver + + ! flow variables, parameter arrays + real(WP), dimension(:,:,:,:), pointer :: Uc, dU !< state variables + real(WP), dimension(:,:,:,:), pointer :: params !< params for evals/flux + + ! CFL numbers + real(WP) :: CFL_x, CFL_y, CFL_z !< global CFL numbers (over dt) + logical :: have_CFL_x, have_CFL_y, have_CFL_z + + ! monitoring quantities + real(WP), dimension(:), pointer :: Umin, Umax !< state variable range + real(WP), dimension(:), pointer :: Uint !< integral of state vars over domain + logical :: have_Urange + + contains + + ! standard interface + procedure :: print => rusanov_print !< output solver to the screen + procedure :: add_bcond !< add a boundary condition + procedure :: get_bcond !< get a boundary condition + procedure :: apply_bcond !< apply all boundary conditions + procedure :: get_cfl !< get maximum CFL + procedure :: get_range !< calculate min/max field values + procedure :: get_min => get_range !< compatibility + procedure :: get_max => get_range !< compatibility + procedure :: calc_dU_x, calc_dU_y, calc_dU_z !< take step + + end type rusanov + + !> declare constructors + interface rusanov; procedure rusanov_from_args; end interface + +contains + + !! standard interface class functions + + !> default constructor for rusanov-type flow solver + function rusanov_from_args(cfg, name, N, P, evals_x, evals_y, evals_z, & + flux_x, flux_y, flux_z) result(this) + use messager, only: die + implicit none + type(rusanov) :: this + class(config), target, intent(in) :: cfg + character(len=*), intent(in), optional :: name + integer, intent(in) :: N, P + procedure(eigenvals_ftype), pointer, intent(in) :: evals_x, evals_y, evals_z + procedure(flux_ftype), pointer, intent(in) :: flux_x, flux_y, flux_z + integer :: imino, imaxo, jmino, jmaxo, kmino, kmaxo, minmino, maxmaxo + + ! check number of ghost cells + if (cfg%no .lt. 1) call die("must have at least one ghost cell") + + ! point to pgrid object + this%cfg => cfg + + ! set the name for the solver + if (present(name)) this%name = trim(adjustl(name)) + + ! system information + this%N = N; this%P = P; + this%evals_x => evals_x; this%evals_y => evals_y; this%evals_z => evals_z; + this%flux_x => flux_x; this%flux_y => flux_y; this%flux_z => flux_z; + + ! boundary condition list + this%nbc = 0 + this%first_bc => NULL() + + ! get array sizes + imino = this%cfg%imino_; imaxo = this%cfg%imaxo_; + jmino = this%cfg%jmino_; jmaxo = this%cfg%jmaxo_; + kmino = this%cfg%kmino_; kmaxo = this%cfg%kmaxo_; + minmino = min(imino, jmino, kmino); maxmaxo = max(imaxo, jmaxo, kmaxo) + + ! allocate flow variables and parameter arrays + allocate(this%params(P,imino:imaxo,jmino:jmaxo,kmino:kmaxo)) + allocate(this%Uc(N,imino:imaxo,jmino:jmaxo,kmino:kmaxo)) + allocate(this%dU(N,imino:imaxo,jmino:jmaxo,kmino:kmaxo)) + + ! allocate velocity masks + allocate(this%vel_mask_x(N), this%vel_mask_y(N), this%vel_mask_z(N)) + + ! initialize CFL computations + this%have_CFL_x = .false. + this%have_CFL_y = .false. + this%have_CFL_z = .false. + + ! initialize monitoring + allocate(this%Umin(N), this%Umax(N), this%Uint(N)) + this%have_Urange = .false. + + end function + + !> solver print function + subroutine rusanov_print(this) + use, intrinsic :: iso_fortran_env, only: output_unit + implicit none + class(rusanov), intent(in) :: this + + if (this%cfg%amRoot) then + write(output_unit,'("rusanov-type solver [",a,"] for config [",a,"]")') & + & trim(this%name), trim(this%cfg%name) + end if + + end subroutine + + !> Add a boundary condition + subroutine add_bcond(this, name, type, locator, dir) + use iterator_class, only: locator_gen_ftype + use string, only: lowercase + use messager, only: die + implicit none + class(rusanov), intent(inout) :: this + character(len=*), intent(in) :: name + integer(1), intent(in) :: type + procedure(locator_gen_ftype) :: locator + character(len=2), intent(in) :: dir + type(rus_bc), pointer :: new_bc + + ! prepare new bcond + allocate(new_bc) + new_bc%name = trim(adjustl(name)) + new_bc%type = type + new_bc%itr = iterator(pg=this%cfg, name=new_bc%name, locator=locator) + select case (lowercase(dir)) + case ('c'); new_bc%dir = 0_1 + case ('-x', 'x-', 'xl'); new_bc%dir = 2_1 + case ('+x', 'x+', 'xr'); new_bc%dir = 1_1 + case ('-y', 'y-', 'yl'); new_bc%dir = 4_1 + case ('+y', 'y+', 'yr'); new_bc%dir = 3_1 + case ('-z', 'z-', 'zl'); new_bc%dir = 6_1 + case ('+z', 'z+', 'zr'); new_bc%dir = 5_1 + case default; call die('[rusanov add_bcond] unknown bcond direction') + end select + + ! insert it in front + new_bc%next => this%first_bc + this%first_bc => new_bc + this%nbc = this%nbc + 1 + + end subroutine add_bcond + + !> get a boundary condition + subroutine get_bcond(this, name, my_bc) + use messager, only: die + implicit none + class(rusanov), intent(inout) :: this + character(len=*), intent(in) :: name + type(rus_bc), pointer, intent(out) :: my_bc + + my_bc => this%first_bc + do while (associated(my_bc)) + if (trim(my_bc%name).eq.trim(name)) return + my_bc => my_bc%next + end do + + call die('[rusanov get_bcond] Boundary condition was not found') + + end subroutine get_bcond + + !> enforce boundary condition + !> note that this interpolates (ensuring ghost cells have the right value), + !> but all wave-related constraints are imposed during the step + subroutine apply_bcond(this, t, dt) + use messager, only: die + implicit none + class(rusanov), intent(inout) :: this + real(WP), intent(in) :: t, dt + logical, dimension(this%N) :: vel_mask + integer :: i, j, k, m, n, iref, jref, kref + type(rus_bc), pointer :: my_bc + + ! Traverse bcond list + my_bc => this%first_bc + do while (associated(my_bc)) + + ! Only processes inside the bcond work here + if (my_bc%itr%amIn) then + + ! Select appropriate action based on the bcond type + select case (my_bc%type) + case (0_1) !< do nothing + case (1_1) !< extrapolate + do m = 1, my_bc%itr%no_ + i = my_bc%itr%map(1,m) + j = my_bc%itr%map(2,m) + k = my_bc%itr%map(3,m) + iref = min(this%cfg%imax, max(this%cfg%imin, i)) + jref = min(this%cfg%jmax, max(this%cfg%jmin, j)) + kref = min(this%cfg%kmax, max(this%cfg%kmin, k)) + this%Uc(:,i,j,k) = this%Uc(:,iref,jref,kref) + end do + case (2_1) !< reflect + do m = 1, my_bc%itr%no_ + i = my_bc%itr%map(1,m) + j = my_bc%itr%map(2,m) + k = my_bc%itr%map(3,m) + iref = min(this%cfg%imax, max(this%cfg%imin, i)) + jref = min(this%cfg%jmax, max(this%cfg%jmin, j)) + kref = min(this%cfg%kmax, max(this%cfg%kmin, k)) + if (i .ne. iref) then + if (i .lt. iref) then + iref = i + 2 * (iref - i) - 1 + else + iref = i - 2 * (i - iref) + 1 + end if + vel_mask(:) = this%vel_mask_x(:) + end if + if (j .ne. jref) then + if (j .lt. jref) then + jref = j + 2 * (jref - j) - 1 + else + jref = j - 2 * (j - jref) + 1 + end if + vel_mask(:) = this%vel_mask_y(:) + end if + if (k .ne. kref) then + if (k .lt. kref) then + kref = k + 2 * (kref - k) - 1 + else + kref = k - 2 * (k - kref) + 1 + end if + vel_mask(:) = this%vel_mask_z(:) + end if + do n = 1, this%N + if (vel_mask(n)) then + this%Uc(n,i,j,k) = - this%Uc(n,iref,jref,kref) + else + this%Uc(n,i,j,k) = + this%Uc(n,iref,jref,kref) + end if + end do + end do + case default + call die('[rusanov apply_bcond] Unknown bcond type') + end select + + ! zero normal velocity? + !masked_type = iand(my_bc%type, 4_1) + !if (masked_type .eq. 4_1) then + ! select case (my_bc%dir) + ! case (0_1) + ! vel_mask(:) = this%vel_mask_x(:) .and. this%vel_mask_y(:) .and. & + ! & this%vel_mask_z(:) + ! case (1_1, 2_1) + ! vel_mask(:) = this%vel_mask_x(:) + ! case (3_1, 4_1) + ! vel_mask(:) = this%vel_mask_y(:) + ! case (5_1, 6_1) + ! vel_mask(:) = this%vel_mask_z(:) + ! end select + ! do m = 1, my_bc%itr%n_ + ! i = my_bc%itr%map(1,m) + ! j = my_bc%itr%map(2,m) + ! k = my_bc%itr%map(3,m) + ! do n = 1, this%N + ! if (vel_mask(n)) this%Uc(n,i,j,k) = 0.0_WP + ! end do + ! end do + !end if + + end if + + ! Sync full fields after each bcond - this should be optimized + call this%cfg%sync(this%Uc) + call this%cfg%sync(this%dU) + + ! Move on to the next bcond + my_bc => my_bc%next + + end do + + end subroutine apply_bcond + + !> get CFL number + subroutine get_cfl(this, dt, cfl) + use messager, only: die + use mpi_f08, only: MPI_ALLREDUCE, MPI_MAX + use parallel, only: MPI_REAL_WP + implicit none + class(rusanov), intent(inout) :: this + real(WP), intent(in) :: dt + real(WP), intent(out) :: cfl + real(WP), dimension(this%N) :: evals + real(WP) :: task_CFL_x, task_CFL_y, task_CFL_z + logical :: have_CFL + integer :: i, j, k, ierrx, ierry, ierrz + + have_CFL = this%have_CFL_x .and. this%have_CFL_y .and. this%have_CFL_z + + if (.not. have_CFL) then + + call this%cfg%sync(this%Uc) + + task_CFL_x = 0.0_WP; task_CFL_y = 0.0_WP; task_CFL_z = 0.0_WP; + do k = this%cfg%kmin_, this%cfg%kmax_ + do j = this%cfg%jmin_, this%cfg%jmax_ + do i = this%cfg%imin_, this%cfg%imax_ + call this%evals_x(this%P, this%N, this%params(:,i,j,k), this%Uc(:,i,j,k), evals) + task_CFL_x = max(task_CFL_x, this%cfg%dxi(i) * maxval(abs(evals))) + call this%evals_y(this%P, this%N, this%params(:,i,j,k), this%Uc(:,i,j,k), evals) + task_CFL_y = max(task_CFL_y, this%cfg%dyi(i) * maxval(abs(evals))) + call this%evals_z(this%P, this%N, this%params(:,i,j,k), this%Uc(:,i,j,k), evals) + task_CFL_z = max(task_CFL_z, this%cfg%dzi(i) * maxval(abs(evals))) + end do + end do + end do + + call MPI_ALLREDUCE(task_CFL_x, this%CFL_x, 1, MPI_REAL_WP, MPI_MAX, & + & this%cfg%comm, ierrx) + call MPI_ALLREDUCE(task_CFL_y, this%CFL_y, 1, MPI_REAL_WP, MPI_MAX, & + & this%cfg%comm, ierry) + call MPI_ALLREDUCE(task_CFL_z, this%CFL_z, 1, MPI_REAL_WP, MPI_MAX, & + & this%cfg%comm, ierrz) + + if (ierrx .ne. 0) call die("error reducing CFL") + if (ierry .ne. 0) call die("error reducing CFL") + if (ierrz .ne. 0) call die("error reducing CFL") + + this%have_CFL_x = .true. + this%have_CFL_y = .true. + this%have_CFL_z = .true. + + end if + + cfl = dt * max(this%CFL_x, this%CFL_y, this%CFL_z) + + end subroutine get_cfl + + !> get min/max field values + subroutine get_range(this) + use messager, only: die + use mpi_f08, only: MPI_ALLREDUCE, MPI_MIN, MPI_MAX, MPI_SUM + use parallel, only: MPI_REAL_WP + implicit none + class(rusanov), intent(inout) :: this + real(WP), dimension(this%N) :: task_Umin, task_Umax, task_Uint + integer :: i, j, k, ierr + + if (.not. this%have_Urange) then + + task_Umin(:) = + huge(1.0_WP) + task_Umax(:) = - huge(1.0_WP) + task_Uint(:) = 0.0_WP + + do k=this%cfg%kmin_,this%cfg%kmax_ + do j=this%cfg%jmin_,this%cfg%jmax_ + do i=this%cfg%imin_,this%cfg%imax_ + task_Umin = min(task_Umin, this%Uc(:,i,j,k)) + task_Umax = max(task_Umax, this%Uc(:,i,j,k)) + task_Uint = task_Uint + this%Uc(:,i,j,k) * this%cfg%dx(i) * & + & this%cfg%dy(j) * this%cfg%dz(k) + end do + end do + end do + + call MPI_ALLREDUCE(task_Umin, this%Umin, this%N, MPI_REAL_WP, MPI_MIN, & + & this%cfg%comm, ierr) + if (ierr .ne. 0) call die("error reducing Umin") + call MPI_ALLREDUCE(task_Umax, this%Umax, this%N, MPI_REAL_WP, MPI_MAX, & + & this%cfg%comm, ierr) + if (ierr .ne. 0) call die("error reducing Umax") + call MPI_ALLREDUCE(task_Uint, this%Uint, this%N, MPI_REAL_WP, MPI_SUM, & + & this%cfg%comm, ierr) + if (ierr .ne. 0) call die("error reducing Uint") + + this%have_Urange = .true. + + end if + + end subroutine get_range + + !> compute dU in x direction + subroutine calc_dU_x(this, dt) + implicit none + class(rusanov), intent(inout) :: this + real(WP), intent(in) :: dt + integer :: j, k + + call this%cfg%sync(this%Uc) + + do k = this%cfg%kmin_, this%cfg%kmax_ + do j = this%cfg%jmin_, this%cfg%jmax_ + call rusanov_1d(this%N, this%P, size(this%Uc, 2), this%flux_x, & + this%evals_x, this%cfg%dx, dt, this%params(:,:,j,k), & + this%Uc(:,:,j,k), this%dU(:,:,j,k)) + end do + end do + + call this%cfg%sync(this%dU) + + this%have_CFL_x = .false.; this%have_Urange = .false.; + + end subroutine calc_dU_x + + !> compute dU in y direction + subroutine calc_dU_y(this, dt) + implicit none + class(rusanov), intent(inout) :: this + real(WP), intent(in) :: dt + integer :: i, k, jmino_, jmaxo_ + real(WP), dimension(:,:), allocatable :: U1d, dU1d, p1d + + if (this%cfg%ny .lt. 2) return + + call this%cfg%sync(this%Uc) + + jmino_ = this%cfg%jmino_; jmaxo_ = this%cfg%jmaxo_; + + allocate(U1d(this%N,jmino_:jmaxo_), p1d(this%P,jmino_:jmaxo_), dU1d(this%N,jmino_:jmaxo_)) + + do k = this%cfg%kmin_, this%cfg%kmax_ + do i = this%cfg%imin_, this%cfg%imax_ + U1d(:,:) = this%Uc(:,i,jmino_:jmaxo_,k) + dU1d(:,:) = this%dU(:,i,jmino_:jmaxo_,k) + p1d(:,:) = this%params(:,i,jmino_:jmaxo_,k) + call rusanov_1d(this%N, this%P, size(U1d, 2), this%flux_y, & + this%evals_y, this%cfg%dy, dt, p1d, U1d, dU1d) + this%dU(:,i,jmino_:jmaxo_,k) = dU1d + end do + end do + + deallocate(U1d, dU1d, p1d) + + call this%cfg%sync(this%dU) + + this%have_CFL_y = .false.; this%have_Urange = .false.; + + end subroutine calc_dU_y + + !> compute dU in z direction + subroutine calc_dU_z(this, dt) + implicit none + class(rusanov), intent(inout) :: this + real(WP), intent(in) :: dt + integer :: i, j, kmino_, kmaxo_ + real(WP), dimension(:,:), allocatable :: U1d, dU1d, p1d + + if (this%cfg%nz .lt. 2) return + + call this%cfg%sync(this%Uc) + + kmino_ = this%cfg%kmino_; kmaxo_ = this%cfg%kmaxo_; + + allocate(U1d(this%N,kmino_:kmaxo_), p1d(this%P,kmino_:kmaxo_), dU1d(this%N,kmino_:kmaxo_)) + + do j = this%cfg%jmin_, this%cfg%jmax_ + do i = this%cfg%imin_, this%cfg%imax_ + U1d(:,:) = this%Uc(:,i,j,kmino_:kmaxo_) + dU1d(:,:) = this%dU(:,i,j,kmino_:kmaxo_) + p1d(:,:) = this%params(:,i,j,kmino_:kmaxo_) + call rusanov_1d(this%N, this%P, size(U1d, 2), this%flux_z, & + this%evals_z, this%cfg%dz, dt, p1d, U1d, dU1d) + this%dU(:,i,j,kmino_:kmaxo_) = dU1d + end do + end do + + deallocate(U1d, dU1d, p1d) + + call this%cfg%sync(this%dU) + + this%have_CFL_z = .false.; this%have_Urange = .false.; + + end subroutine calc_dU_z + + ! requires one ghost cell + ! left as an independent function, but called is by class + pure subroutine rusanov_1d(N, P, M, flux, evals, dxs, dt, params, U, dU) + implicit none + integer, intent(in) :: N, P, M + procedure(eigenvals_ftype), pointer, intent(in) :: evals + procedure(flux_ftype), pointer, intent(in) :: flux + real(WP), dimension(M), intent(in) :: dxs ! size M + real(WP), intent(in) :: dt + real(WP), dimension(P,M), intent(in) :: params ! size (P, M) + real(WP), dimension(N,M), intent(in) :: U ! size (N, M) + real(WP), dimension(N,M), intent(inout) :: dU ! size (N, M) + real(WP), dimension(N) :: F, Fl, Fr, nudiff, laml, lamr + integer :: j + + dU(:,:) = 0.0_WP + + call flux(P, N, params(:,1), U(:,1), Fl) + call evals(P, N, params(:,1), U(:,1), laml) + laml(:) = abs(laml(:)) + + do j = 2, M + + call flux(P, N, params(:,j), U(:,j), Fr) + call evals(P, N, params(:,j), U(:,j), lamr) + lamr(:) = abs(lamr(:)) + + nudiff = (U(:,j) - U(:,j-1)) * max(maxval(laml), maxval(lamr)) + + ! not the right division by dx, but whatever + F = dt * 0.5_WP / (0.5_WP * (dxs(j-1) + dxs(j))) * (Fr + Fl - nudiff) + + dU(:,j-1) = dU(:,j-1) - F + dU(:,j ) = dU(:,j ) + F + + Fl(:) = Fr(:); laml(:) = lamr(:); + + end do + + end subroutine rusanov_1d + +end module rusanov_class +