output_partial_domain for post-processing (#768)

Co-authored-by: Spencer Bryngelson <[email protected]>

docs/documentation/case.md

@@ -522,6 +522,9 @@ To restart the simulation from $k$-th time step, see [Restarting Cases](running.
 | `probe_wrt`          | Logical | Write the flow chosen probes data files for each time step	|
 | `num_probes`         | Integer | Number of probes	|
 | `probe(i)%[x,y,z]`   | Real	 | Coordinates of probe $i$	|
+| `output_partial_domain` | Logical | Output part of the domain |
+| `[x,y,z]_output%beg` | Real    | Beginning of the output domain in the [x,y,z]-direction |
+| `[x,y,z]_output%end` | Real    | End of the output domain in the [x,y,z]-direction |
 
 The table lists formatted database output parameters. The parameters define variables that are outputted from simulation and file types and formats of data as well as options for post-processing.
 
@@ -549,6 +552,11 @@ If `file_per_process` is true, then pre_process, simulation, and post_process mu
 
 - `probe_wrt` activates the output of state variables at coordinates specified by `probe(i)%[x;y,z]`.
 
+- `output_partial_domain` activates the output of part of the domain specified by `[x,y,z]_output%beg` and `[x,y,z]_output%end`.
+This is useful for large domains where only a portion of the domain is of interest.
+It is not supported when `precision = 1` and `format = 1`. 
+It also cannot be enabled with `flux_wrt`, `heat_ratio_wrt`, `pres_inf_wrt`, `c_wrt`, `omega_wrt`, `ib`, `schlieren_wrt`, or `qm_wrt`.
+
 ### 8. Acoustic Source {#acoustic-source}
 
 | Parameter                             | Type    | Description |

src/post_process/m_checker.fpp

@@ -26,6 +26,7 @@ contains
     subroutine s_check_inputs
 
         call s_check_inputs_output_format
+        call s_check_inputs_partial_domain
         call s_check_inputs_partial_density
         call s_check_inputs_velocity
         call s_check_inputs_flux_limiter
@@ -43,6 +44,21 @@ contains
         @:PROHIBIT(precision /= 1 .and. precision /= 2)
     end subroutine s_check_inputs_output_format
 
+    !> Checks constraints on partial domain parameters
+    subroutine s_check_inputs_partial_domain
+        @:PROHIBIT(output_partial_domain .and. format == 1)
+        @:PROHIBIT(output_partial_domain .and. precision == 1)
+        @:PROHIBIT(output_partial_domain .and. any([flux_wrt, heat_ratio_wrt, pres_inf_wrt, c_wrt, schlieren_wrt, qm_wrt, ib, any(omega_wrt)]))
+
+        @:PROHIBIT(output_partial_domain .and. (f_is_default(x_output%beg) .or. f_is_default(x_output%end)))
+        @:PROHIBIT(output_partial_domain .and. n /= 0 .and. (f_is_default(y_output%beg) .or. f_is_default(y_output%end)))
+        @:PROHIBIT(output_partial_domain .and. p /= 0 .and. (f_is_default(z_output%beg) .or. f_is_default(z_output%end)))
+
+        #:for X in ['x', 'y', 'z']
+            @:PROHIBIT(${X}$_output%beg > ${X}$_output%end)
+        #:endfor
+    end subroutine s_check_inputs_partial_domain
+
     !> Checks constraints on partial density parameters
     subroutine s_check_inputs_partial_density
         character(len=5) :: iStr

src/post_process/m_data_output.fpp

@@ -28,6 +28,7 @@ module m_data_output
     implicit none
 
     private; public :: s_initialize_data_output_module, &
+ s_define_output_region, &
  s_open_formatted_database_file, &
  s_open_intf_data_file, &
  s_open_energy_data_file, &
@@ -365,6 +366,9 @@ contains
             ! Pressure
             if (pres_wrt .or. prim_vars_wrt) dbvars = dbvars + 1
 
+            ! Elastic stresses
+            if (hypoelasticity) dbvars = dbvars + (num_dims*(num_dims + 1))/2
+
             ! Volume fraction(s)
             if ((model_eqns == 2) .or. (model_eqns == 3)) then
 
@@ -430,6 +434,42 @@ contains
 
     end subroutine s_initialize_data_output_module ! --------------------------
 
+    subroutine s_define_output_region
+
+        integer :: i
+        integer :: lower_bound, upper_bound
+
+        #:for X, M in [('x', 'm'), ('y', 'n'), ('z', 'p')]
+
+            if (${M}$ == 0) return ! Early return for y or z if simulation is 1D or 2D
+
+            lower_bound = -offset_${X}$%beg
+            upper_bound = ${M}$+offset_${X}$%end
+
+            do i = lower_bound, upper_bound
+                if (${X}$_cc(i) > ${X}$_output%beg) then
+                    ${X}$_output_idx%beg = i + offset_${X}$%beg
+                    exit
+                end if
+            end do
+
+            do i = upper_bound, lower_bound, -1
+                if (${X}$_cc(i) < ${X}$_output%end) then
+                    ${X}$_output_idx%end = i + offset_${X}$%beg
+                    exit
+                end if
+            end do
+
+            ! If no grid points are within the output region
+            if ((${X}$_cc(lower_bound) > ${X}$_output%end) .or. (${X}$_cc(upper_bound) < ${X}$_output%beg)) then
+                ${X}$_output_idx%beg = 0
+                ${X}$_output_idx%end = 0
+            end if
+
+        #:endfor
+
+    end subroutine s_define_output_region
+
     subroutine s_open_formatted_database_file(t_step) ! --------------------
         ! Description: This subroutine opens a new formatted database file, or
         !              replaces an old one, and readies it for the data storage
@@ -520,7 +560,14 @@ contains
             ! file by describing in it the dimensionality of post-processed
             ! data as well as the total number of flow variable(s) that will
             ! eventually be stored in it
-            write (dbfile) m, n, p, dbvars
+            if (output_partial_domain) then
+                write (dbfile) x_output_idx%end - x_output_idx%beg, &
+                    y_output_idx%end - y_output_idx%beg, &
+                    z_output_idx%end - z_output_idx%beg, &
+                    dbvars
+            else
+                write (dbfile) m, n, p, dbvars
+            end if
 
             ! Next, analogous steps to the ones above are carried out by the
             ! root process to create and setup the formatted database master
@@ -539,7 +586,11 @@ contains
                                      '. Exiting ...')
                 end if
 
-                write (dbroot) m_root, 0, 0, dbvars
+                if (output_partial_domain) then
+                    write (dbroot) x_output_idx%end - x_output_idx%beg, 0, 0, dbvars
+                else
+                    write (dbroot) m_root, 0, 0, dbvars
+                end if
 
             end if
 
@@ -729,15 +780,26 @@ contains
                         real(y_cb, sp), &
                         real(z_cb, sp)
                 else
-                    write (dbfile) x_cb, y_cb, z_cb
+                    if (output_partial_domain) then
+                        write (dbfile) x_cb(x_output_idx%beg - 1:x_output_idx%end), &
+                            y_cb(y_output_idx%beg - 1:y_output_idx%end), &
+                            z_cb(z_output_idx%beg - 1:z_output_idx%end)
+                    else
+                        write (dbfile) x_cb, y_cb, z_cb
+                    end if
                 end if
 
             elseif (n > 0) then
                 if (precision == 1) then
                     write (dbfile) real(x_cb, sp), &
                         real(y_cb, sp)
                 else
-                    write (dbfile) x_cb, y_cb
+                    if (output_partial_domain) then
+                        write (dbfile) x_cb(x_output_idx%beg - 1:x_output_idx%end), &
+                            y_cb(y_output_idx%beg - 1:y_output_idx%end)
+                    else
+                        write (dbfile) x_cb, y_cb
+                    end if
                 end if
 
                 ! One-dimensional local grid data is written to the formatted
@@ -746,9 +808,13 @@ contains
             else
 
                 if (precision == 1) then
-                    write (dbfile) real(x_cb, wp)
+                    write (dbfile) real(x_cb, sp)
                 else
-                    write (dbfile) x_cb
+                    if (output_partial_domain) then
+                        write (dbfile) x_cb(x_output_idx%beg - 1:x_output_idx%end)
+                    else
+                        write (dbfile) x_cb
+                    end if
                 end if
 
                 if (num_procs > 1) then
@@ -761,7 +827,11 @@ contains
                     if (precision == 1) then
                         write (dbroot) real(x_root_cb, wp)
                     else
-                        write (dbroot) x_root_cb
+                        if (output_partial_domain) then
+                            write (dbroot) x_root_cb(x_output_idx%beg - 1:x_output_idx%end)
+                        else
+                            write (dbroot) x_root_cb
+                        end if
                     end if
                 end if
 

src/post_process/m_global_parameters.fpp

@@ -197,6 +197,11 @@ module m_global_parameters
 
     integer :: precision !< Floating point precision of the database file(s)
 
+    logical :: output_partial_domain !< Specify portion of domain to output for post-processing
+
+    type(bounds_info) :: x_output, y_output, z_output !< Portion of domain to output for post-processing
+    type(int_bounds_info) :: x_output_idx, y_output_idx, z_output_idx !< Indices of domain to output for post-processing
+
     !> @name Size of the ghost zone layer in the x-, y- and z-coordinate directions.
     !! The definition of the ghost zone layers is only necessary when using the
     !! Silo database file format in multidimensions. These zones provide VisIt
@@ -432,6 +437,15 @@ contains
         ! IBM
         num_ibs = dflt_int
 
+        ! Output partial domain
+        output_partial_domain = .false.
+        x_output%beg = dflt_real
+        x_output%end = dflt_real
+        y_output%beg = dflt_real
+        y_output%end = dflt_real
+        z_output%beg = dflt_real
+        z_output%end = dflt_real
+
     end subroutine s_assign_default_values_to_user_inputs
 
     !>  Computation of parameters, allocation procedures, and/or
@@ -692,6 +706,15 @@ contains
             species_idx%end = 1
         end if
 
+        if (output_partial_domain) then
+            x_output_idx%beg = 0
+            x_output_idx%end = 0
+            y_output_idx%beg = 0
+            y_output_idx%end = 0
+            z_output_idx%beg = 0
+            z_output_idx%end = 0
+        end if
+
         momxb = mom_idx%beg
         momxe = mom_idx%end
         advxb = adv_idx%beg

src/post_process/m_mpi_proxy.fpp

@@ -170,7 +170,7 @@ contains
             & 'prim_vars_wrt', 'c_wrt', 'qm_wrt','schlieren_wrt', 'bubbles_euler', 'qbmm',   &
             & 'polytropic', 'polydisperse', 'file_per_process', 'relax', 'cf_wrt',     &
             & 'adv_n', 'ib', 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt',          &
-            & 'surface_tension', 'hyperelasticity', 'bubbles_lagrange', 'rkck_adap_dt']
+            & 'surface_tension', 'hyperelasticity', 'bubbles_lagrange', 'rkck_adap_dt', 'output_partial_domain']
             call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 
@@ -191,7 +191,9 @@ contains
         end do
 
         #:for VAR in [ 'pref', 'rhoref', 'R0ref', 'poly_sigma', 'Web', 'Ca', &
-            & 'Re_inv', 'sigma', 't_save', 't_stop' ]
+            & 'Re_inv', 'sigma', 't_save', 't_stop', &
+            & 'x_output%beg', 'x_output%end', 'y_output%beg', &
+            & 'y_output%end', 'z_output%beg', 'z_output%end' ]
             call MPI_BCAST(${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
         #:endfor
         call MPI_BCAST(schlieren_alpha(1), num_fluids_max, mpi_p, 0, MPI_COMM_WORLD, ierr)