From f45888f9933aac983788f76a23e3dfd55fb92e27 Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Thu, 7 Nov 2024 11:33:00 -0500 Subject: [PATCH 1/9] Add viscous and surface_tension logicals (#688) Co-authored-by: Spencer Bryngelson --- benchmarks/viscous_weno5_sgb_acoustic/case.py | 23 ++- docs/documentation/case.md | 74 +++---- examples/2D_TaylorGreenVortex/case.py | 3 +- examples/2D_bubbly_steady_shock/case.py | 1 - examples/2D_ibm/case.py | 1 + examples/2D_ibm_cfl_dt/case.py | 1 + examples/2D_ibm_multiphase/case.py | 3 - examples/2D_laplace_pressure_jump/case.py | 12 +- examples/2D_lid_driven_cavity/case.py | 3 +- examples/2D_mixing_artificial_Ma/case.py | 1 + examples/2D_rayleigh_taylor/case.py | 1 + examples/2D_shearlayer/case.py | 1 + examples/2D_viscous/case.py | 1 + examples/3D_TaylorGreenVortex/case.py | 1 + examples/3D_ibm_bowshock/case.py | 3 +- examples/3D_rayleigh_taylor/case.py | 1 + examples/3D_recovering_sphere/case.py | 3 +- examples/3D_turb_mixing/case.py | 1 + src/common/m_checker_common.fpp | 11 +- src/common/m_mpi_common.fpp | 4 +- src/common/m_variables_conversion.fpp | 12 +- src/post_process/m_checker.fpp | 2 +- src/post_process/m_global_parameters.fpp | 6 +- src/post_process/m_start_up.f90 | 2 +- src/pre_process/m_assign_variables.fpp | 2 +- src/pre_process/m_global_parameters.fpp | 6 +- src/pre_process/m_mpi_proxy.fpp | 2 +- src/pre_process/m_start_up.fpp | 3 +- src/simulation/m_checker.fpp | 7 + src/simulation/m_data_output.fpp | 30 +-- src/simulation/m_global_parameters.fpp | 27 ++- src/simulation/m_mpi_proxy.fpp | 7 +- src/simulation/m_rhs.fpp | 190 +++++++++--------- src/simulation/m_riemann_solvers.fpp | 112 +++++------ src/simulation/m_sim_helpers.f90 | 12 +- src/simulation/m_start_up.fpp | 13 +- src/simulation/m_time_steppers.fpp | 2 +- src/simulation/m_viscous.fpp | 20 +- toolchain/mfc/run/case_dicts.py | 6 +- toolchain/mfc/test/cases.py | 18 +- 40 files changed, 331 insertions(+), 297 deletions(-) diff --git a/benchmarks/viscous_weno5_sgb_acoustic/case.py b/benchmarks/viscous_weno5_sgb_acoustic/case.py index 4017fa955..ccdde880c 100644 --- a/benchmarks/viscous_weno5_sgb_acoustic/case.py +++ b/benchmarks/viscous_weno5_sgb_acoustic/case.py @@ -82,7 +82,7 @@ cact = 1475. t0 = x0/c0 -nbubbles = 1 +nbubbles = 1 myr0 = R0ref cfl = 0.01 @@ -98,7 +98,7 @@ # Logistics ================================================ 'run_time_info' : 'F', # ========================================================== - + # Computational Domain Parameters ========================== 'x_domain%beg' : -10.E-03/x0, 'x_domain%end' : 10.E-03/x0, @@ -127,7 +127,7 @@ 'time_stepper' : 3, 'weno_order' : 5, 'weno_eps' : 1.E-16, - 'weno_Re_flux' : 'F', + 'weno_Re_flux' : 'F', 'weno_avg' : 'F', 'mapped_weno' : 'T', 'null_weights' : 'F', @@ -141,15 +141,16 @@ 'bc_y%end' : -3, 'bc_z%beg' : -3, 'bc_z%end' : -3, + 'viscous' : 'T', # ========================================================== - + # Formatted Database Files Structure Parameters ============ 'format' : 1, 'precision' : 2, 'prim_vars_wrt' :'T', 'parallel_io' :'T', - # ========================================================== - + # ========================================================== + # Patch 1 _ Background ===================================== 'patch_icpp(1)%geometry' : 9, 'patch_icpp(1)%x_centroid' : 0., @@ -167,7 +168,7 @@ 'patch_icpp(1)%r0' : 1., 'patch_icpp(1)%v0' : 0.0E+00, # ========================================================== - + # Patch 2 Screen =========================================== 'patch_icpp(2)%geometry' : 9, 'patch_icpp(2)%x_centroid' : 0., @@ -186,7 +187,7 @@ 'patch_icpp(2)%r0' : 1., 'patch_icpp(2)%v0' : 0.0E+00, # ========================================================== - + # Fluids Physical Parameters =============================== # Surrounding liquid 'fluid_pp(1)%gamma' : 1.E+00/(n_tait-1.E+00), @@ -208,12 +209,12 @@ 'fluid_pp(2)%mu_v' : mu_n, 'fluid_pp(2)%k_v' : k_n, # ========================================================== - + # Non-polytropic gas compression model AND/OR Tait EOS ===== 'pref' : p0, 'rhoref' : rho0, # ========================================================== - + # Bubbles ================================================== 'bubbles' : 'T', 'bubble_model' : 3, @@ -227,7 +228,7 @@ 'Web' : We, 'Re_inv' : Re_inv, # ========================================================== - + # Acoustic source ========================================== 'acoustic_source' : 'T', 'num_source' : 1, diff --git a/docs/documentation/case.md b/docs/documentation/case.md index d2c247a9a..7903cb2a1 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -371,6 +371,8 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r | `n_start` | Integer | Save file from which to start simulation | | `t_save` | Real | Time duration between data output | | `t_stop` | Real | Simulation stop time | +| `surface_tension` | Logical | Activate surface tension | +| `viscous` | Logical | Activate viscosity | - \* Options that work only with `model_eqns = 2`. - † Options that work only with ``cyl_coord = 'F'``. @@ -445,33 +447,37 @@ If this option is false, velocity gradient is computed using finite difference s - `weno_avg` it activates the arithmetic average of the left and right, WENO-reconstructed, cell-boundary values. This option requires `weno_Re_flux` to be true because cell boundary values are only utilized when employing the scalar divergence method in the computation of velocity gradients. +- `surface_tension` activates surface tension when set to ``'T'``. Requires `sigma` to be set. + +- `viscous` activates viscosity when set to ``'T'``. Requires `Re(1)` and `Re(2)` to be set. + #### Constant Time-Stepping -- `dt` specifies the constant time step size that is used in simulation. -The value of `dt` needs to be sufficiently small such that the Courant-Friedrichs-Lewy (CFL) condition is satisfied. +- `dt` specifies the constant time step size used in the simulation. +The value of `dt` needs to be sufficiently small to satisfy the Courant-Friedrichs-Lewy (CFL) condition. -- `t_step_start` and `t_step_end` define the time steps at which simulation starts and ends, respectively. +- `t_step_start` and `t_step_end` define the time steps at which the simulation starts and ends. `t_step_save` is the time step interval for data output during simulation. To newly start the simulation, set `t_step_start = 0`. -To restart simulation from $k$-th time step, set `t_step_start = k`, see [Restarting Cases](running.md#restarting-cases). +To restart the simulation from $k$-th time step, set `t_step_start = k`; see [Restarting Cases](running.md#restarting-cases). ##### Adaptive Time-Stepping - `cfl_adap_dt` enables adaptive time stepping with a constant CFL when true -- `cfl_const_dt` enables constant dt time-stepping where dt results in a specified CFL for the initial condition +- `cfl_const_dt` enables constant `dt` time-stepping where `dt` results in a specified CFL for the initial condition - `cfl_target` specifies the target CFL value - `n_start` specifies the save file to start at -- `t_save` specifies the time interval between data output during simulation +- `t_save` specifies the time interval between data output during the simulation - `t_stop` specifies at what time the simulation should stop To newly start the simulation, set `n_start = 0`. -To restart simulation from $k$-th time step, see [Restarting Cases](running.md#restarting-cases). +To restart the simulation from $k$-th time step, see [Restarting Cases](running.md#restarting-cases). ### 7. Formatted Output @@ -512,23 +518,23 @@ The table lists formatted database output parameters. The parameters define vari - `parallel_io` activates parallel input/output (I/O) of data files. It is highly recommended to activate this option in a parallel environment. With parallel I/O, MFC inputs and outputs a single file throughout pre-process, simulation, and post-process, regardless of the number of processors used. -Parallel I/O enables the use of different number of processors in each of the processes (i.e., simulation data generated using 1000 processors can be post-processed using a single processor). +Parallel I/O enables the use of different numbers of processors in each process (e.g., simulation data generated using 1000 processors can be post-processed using a single processor). - `file_per_process` deactivates shared file MPI-IO and activates file per process MPI-IO. The default behavior is to use a shared file. -File per process is useful when running on 10's of thousands of ranks. +File per process is useful when running on >10K ranks. If `file_per_process` is true, then pre_process, simulation, and post_process must be run with the same number of ranks. -- `cons_vars_wrt` and `prim_vars_wrt` activate output of conservative and primitive state variables into the database, respectively. +- `cons_vars_wrt` and `prim_vars_wrt` activate the output of conservative and primitive state variables into the database. -- `[variable's name]_wrt` activates output of the each specified variable into the database. +- `[variable's name]_wrt` activates the output of each specified variable into the database. - `schlieren_alpha(i)` specifies the intensity of the numerical Schlieren of $i$-th component. -- `fd_order` specifies the order of the finite difference scheme that is used to compute the vorticity from the velocity field and the numerical schlieren from the density field by an integer of 1, 2, and 4. -`fd_order = 1`, `2`, and `4` correspond to the first, second, and fourth-order finite difference schemes, respectively. +- `fd_order` specifies the order of the finite difference scheme used to compute the vorticity from the velocity field and the numerical schlieren from the density field using an integer of 1, 2, and 4. +`fd_order = 1`, `2`, and `4` correspond to the first, second, and fourth-order finite difference schemes. -- `probe_wrt` activates output of state variables at coordinates specified by `probe(i)%[x;y,z]`. +- `probe_wrt` activates the output of state variables at coordinates specified by `probe(i)%[x;y,z]`. ### 8. Acoustic Source {#acoustic-source} @@ -590,15 +596,15 @@ Details of the transducer acoustic source model can be found in [Maeda and Colon - `%%foc_length` specifies the focal length of the transducer for transducer waves. It is the distance from the transducer to the focal point. -- `%%aperture` specifies the aperture of the transducer. It is the diameter of the projection of the transducer arc onto the y-axis (2D) or spherical cap onto the y-z plane (3D). To simulate a transducer enclosing half of the circle/sphere, set the aperture to double the focal length. For transducer array, it is the total aperture of the array. +- `%%aperture` specifies the aperture of the transducer. It is the diameter of the projection of the transducer arc onto the y-axis (2D) or spherical cap onto the y-z plane (3D). Set the aperture to double the focal length to simulate a transducer enclosing half of the circle/sphere. For the transducer array, it is the total aperture of the array. - `%%num_elements` specifies the number of transducer elements in a transducer array. -- `%%element_on` specifies the element number of the transducer array that is on. The element number starts from 1. If all elements are on, set `%%element_on` to 0. +- `%%element_on` specifies the element number of the transducer array that is on. The element number starts from 1, if all elements are on, set `%%element_on` to 0. -- `%%element_spacing_angle` specifies the spacing angle between adjacent transducer in radians. The total aperture (`%%aperture`) is set, so each transducer element is smaller if `%%element_spacing_angle` is larger. +- `%%element_spacing_angle` specifies the spacing angle between adjacent transducers in radians. The total aperture (`%%aperture`) is set, so each transducer element is smaller if `%%element_spacing_angle` is larger. -- `%%element_polygon_ratio` specifies the ratio of the polygon side length to the aperture diameter of each transducer element in a circular 3D transducer array. The polygon side length is calculated by using the total aperture (`%%aperture`) as the circumcicle diameter, and `%%num_elements` as the number of sides of the polygon. The ratio is used specify the aperture size of each transducer element in the array, as a ratio of the total aperture. +- `%%element_polygon_ratio` specifies the ratio of the polygon side length to the aperture diameter of each transducer element in a circular 3D transducer array. The polygon side length is calculated by using the total aperture (`%%aperture`) as the circumcircle diameter and `%%num_elements` as the number of sides of the polygon. The ratio is used to specify the aperture size of each transducer element in the array as a ratio of the total aperture. - `%%rotate_angle` specifies the rotation angle of the 3D circular transducer array along the x-axis (principal axis). It is optional and defaults to 0. @@ -626,16 +632,16 @@ Details of the transducer acoustic source model can be found in [Maeda and Colon | `mu_v` † | Real | Viscosity | | `k_v` † | Real | Thermal conductivity | | `qbmm` | Logical | Quadrature by method of moments| -| `dist_type` | Integer | Joint probability density function for bubble radius and velocity (only when qbmm is true)| -| `sigR` | Real | Standard deviation for probability density function of bubble radius (only when qbmm is true) | -| `sigV` | Real | Standard deviation for probability density function of bubble velocity (only when qbmm is true) | -| `rhoRV` | Real | Correlation coefficient for joint probability density function of bubble radius and velocity (only when qbmm is true) | +| `dist_type` | Integer | Joint probability density function for bubble radius and velocity (only for ``qbmm = 'T'``)| +| `sigR` | Real | Standard deviation for the probability density function of bubble radius (only for ``qbmm = 'T'``) | +| `sigV` | Real | Standard deviation for the probability density function of bubble velocity (only for ``qbmm = 'T'``) | +| `rhoRV` | Real | Correlation coefficient for the joint probability density function of bubble radius and velocity (only for ``qbmm = 'T'``) | -These options work only for gas-liquid two component flows. +These options work only for gas-liquid two-component flows. Component indexes are required to be 1 for liquid and 2 for gas. -- \* These parameters should be pretended with patch index $1$ that is filled with liquid: `fluid_pp(1)%`. -- † These parameters should be pretended with patch indexes that are respectively filled with liquid and gas: `fluid_pp(1)%` and `fluid_pp(2)%`. +- \* These parameters should be prepended with patch index $1$ that is filled with liquid: `fluid_pp(1)%`. +- † These parameters should be prepended with patch indexes filled with liquid and gas: `fluid_pp(1)%` and `fluid_pp(2)%`. This table lists the ensemble-averaged bubble model parameters. @@ -645,9 +651,9 @@ This table lists the ensemble-averaged bubble model parameters. `bubble_model = 1`, `2`, and `3` correspond to the Gilmore, Keller-Miksis, and Rayleigh-Plesset models. - `polytropic` activates polytropic gas compression in the bubble. -When `polytropic` is set `False`, the gas compression is modeled as non-polytropic due to heat and mass transfer across the bubble wall with constant heat and mass transfer coefficients based on ([Preston et al., 2007](references.md#Preston07)). +When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due to heat and mass transfer across the bubble wall with constant heat and mass transfer coefficients based on ([Preston et al., 2007](references.md#Preston07)). -- `polydisperse` activates polydispersity in the bubble model by means of a probability density function (PDF) of the equiliibrium bubble radius. +- `polydisperse` activates polydispersity in the bubble model through a probability density function (PDF) of the equilibrium bubble radius. - `thermal` specifies a model for heat transfer across the bubble interface by an integer from 1 through 3. `thermal = 1`, `2`, and `3` correspond to no heat transfer (adiabatic gas compression), isothermal heat transfer, and heat transfer with a constant heat transfer coefficient based on [Preston et al., 2007](references.md#Preston07), respectively. @@ -672,11 +678,11 @@ Implementation of the parameters into the model follow [Ando (2010)](references. - `dist_type` specifies the initial joint PDF of initial bubble radius and bubble velocity required in qbmm. `dist_type = 1` and `2` correspond to binormal and lognormal-normal distributions respectively. -- `sigR` specifies the standard deviation of the PDF of bubble radius required in qbmm. +- `sigR` specifies the standard deviation of the PDF of bubble radius required in the QBMM feature. -- `sigV` specifies the standard deviation of the PDF of bubble velocity required in qbmm. +- `sigV` specifies the standard deviation of the PDF of bubble velocity required in the QBMM feature. -- `rhoRV` specifies the correlation coefficient of the joint PDF of bubble radius and bubble velocity required in qbmm. +- `rhoRV` specifies the correlation coefficient of the joint PDF of bubble radius and bubble velocity required in the QBMM feature. ### 10. Velocity Field Setup @@ -705,7 +711,7 @@ The parameters are optionally used to define initial velocity profiles and pertu - `perturb_sph_fluid` specifies the fluid component whose partial density is to be perturbed. -- `mixlayer_vel_profile` activates setting of the mean streamwise velocity to hyperbolic tangent profile. This option works only for 2D and 3D cases. +- `mixlayer_vel_profile` activates setting the mean streamwise velocity to a hyperbolic tangent profile. This option works only for 2D and 3D cases. - `mixlayer_vel_coef` is a parameter for the hyperbolic tangent profile of a mixing layer when `mixlayer_vel_profile = 'T'`. The mean streamwise velocity profile is given as: @@ -727,7 +733,7 @@ $$ u = patch\_icpp(1)\%vel(1) * tanh(y\_cc * mixlayer\_vel\_profile) $$ - `relax_model` Specifies the phase change model to be used: [5] enables pT-equilibrium, while [6] activates pTg-equilibrium (if criteria are met). -- `palpha_eps` Specifies the tolerance used for the Newton Solvers used in the pT-equilibrium model. +- `palpha_eps` Specifies the tolerance for the Newton Solvers used in the pT-equilibrium model. - `ptgalpha_eps` Specifies the tolerance used for the Newton Solvers used in the pTg-equilibrium model. @@ -855,7 +861,7 @@ Each patch requires a different set of parameters, which are also listed in this | 10 | Annular Transducer Array | 2D-Axisym | #9 requirements | | 11 | Circular Transducer Array | 3D | #7 requirements, `%%element_polygon_ratio`, and `%%rotate_angle` | -Details of the required parameters for each acoustic support type are listed in [Acoustic Source](#acoustic-source). +The required parameters for each acoustic support type are listed in [Acoustic Source](#acoustic-source). The acoustic support number (`#`) corresponds to the acoustic support type `Acoustic(i)%%support`, where $i$ is the acoustic source index. For each `%%parameter`, prepend the parameter with `acoustic(i)%`. @@ -864,7 +870,7 @@ Additional requirements for all acoustic support types: - `num_source` must be set to the total number of acoustic sources. -- `%%support` must be set to the acoustic support number listed in the table. +- `%%support` must be set to the acoustic support number in the table. - `%%dipole` is only supported for planar sources. diff --git a/examples/2D_TaylorGreenVortex/case.py b/examples/2D_TaylorGreenVortex/case.py index 606b09ccf..a0b6bbebe 100644 --- a/examples/2D_TaylorGreenVortex/case.py +++ b/examples/2D_TaylorGreenVortex/case.py @@ -53,6 +53,7 @@ 'bc_x%end' : -1, 'bc_y%beg' : -1, 'bc_y%end' : -1, + 'viscous' : 'T', # ========================================================== # Formatted Database Files Structure Parameters ============ @@ -85,4 +86,4 @@ 'fluid_pp(1)%Re(1)' : 1/Mu, # ========================================================== })) -# ============================================================================== \ No newline at end of file +# ============================================================================== diff --git a/examples/2D_bubbly_steady_shock/case.py b/examples/2D_bubbly_steady_shock/case.py index 0e6b93df1..083d40c51 100644 --- a/examples/2D_bubbly_steady_shock/case.py +++ b/examples/2D_bubbly_steady_shock/case.py @@ -192,7 +192,6 @@ 'fluid_pp(1)%M_v' : M_v, 'fluid_pp(1)%mu_v' : mu_v, 'fluid_pp(1)%k_v' : k_v, - #'fluid_pp(1)%Re(1)' : 80000, # Last fluid_pp is always reserved for bubble gas state === # if applicable ========================================== diff --git a/examples/2D_ibm/case.py b/examples/2D_ibm/case.py index e21319e92..289b3633c 100644 --- a/examples/2D_ibm/case.py +++ b/examples/2D_ibm/case.py @@ -62,6 +62,7 @@ # Set IB to True and add 1 patch 'ib' : 'T', 'num_ibs' : 1, + 'viscous' : 'T', # ========================================================================== # Formatted Database Files Structure Parameters ============================ diff --git a/examples/2D_ibm_cfl_dt/case.py b/examples/2D_ibm_cfl_dt/case.py index 199959ce4..5c4326d21 100644 --- a/examples/2D_ibm_cfl_dt/case.py +++ b/examples/2D_ibm_cfl_dt/case.py @@ -64,6 +64,7 @@ # Set IB to True and add 1 patch 'ib' : 'T', 'num_ibs' : 1, + 'viscous' : 'T', # ========================================================================== # Formatted Database Files Structure Parameters ============================ diff --git a/examples/2D_ibm_multiphase/case.py b/examples/2D_ibm_multiphase/case.py index 689d52ea1..c4c457709 100644 --- a/examples/2D_ibm_multiphase/case.py +++ b/examples/2D_ibm_multiphase/case.py @@ -100,10 +100,7 @@ # Specify 2 fluids 'fluid_pp(1)%gamma' : 1.E+00/(gam_a-1.E+00), # 2.50(Not 1.40) 'fluid_pp(1)%pi_inf' : 0, - # 'fluid_pp(1)%Re(1)' : 2500000, - # 'fluid_pp(1)%Re(2)' : 1.0e+6, 'fluid_pp(2)%gamma' : 1.E+00/(gam_b-1.E+00), # 2.50(Not 1.40) 'fluid_pp(2)%pi_inf' : 0, - # 'fluid_pp(2)%Re(1)' : 2500000, # ========================================================================== })) diff --git a/examples/2D_laplace_pressure_jump/case.py b/examples/2D_laplace_pressure_jump/case.py index 508d278e5..97217708d 100644 --- a/examples/2D_laplace_pressure_jump/case.py +++ b/examples/2D_laplace_pressure_jump/case.py @@ -43,8 +43,6 @@ 't_step_start' : 0, 't_step_stop' : 100000, 't_step_save' : 1000, - #'t_step_stop' : 100, - #'t_step_save' : 100, # ======================================= # Simulation Algorithm ================== @@ -53,9 +51,6 @@ 'mixture_err' : 'T', 'mpp_lim' : 'F', 'time_stepper' : 3, - #'recon_type' : 1, - #'muscl_order' : 2, - #'muscl_lim' : 2, 'weno_order' : 5, 'avg_state' : 2, 'weno_eps' : 1e-16, @@ -72,6 +67,7 @@ 'num_patches' : 2, 'num_fluids' : 2, 'weno_avg' : 'T', + 'surface_tension' : 'T', # ======================================= # Database Structure Parameters ========= @@ -84,21 +80,15 @@ # ======================================= 'sigma' : 8, - #'flux_lim' : 2, - #'flux_wrt(1)' : 'T', - #'flux_wrt(2)' : 'T', - #'cf_grad_wrt' : 'T', # Fluid Parameters (Water) ============== 'fluid_pp(1)%gamma' : 1.E+00/(2.1E+00-1.E+00), 'fluid_pp(1)%pi_inf' : 2.1E+00*1.E+06/(2.1E+00-1.E+00), - #'fluid_pp(1)%Re(1)' : 1.e3, # ======================================= # Fluid Parameters (Gas) ================ 'fluid_pp(2)%gamma' : 1.E+00/(1.4E+00-1.E+00), 'fluid_pp(2)%pi_inf' : 0.E+00, - #'fluid_pp(2)%Re(1)' : 1.81e5, # ======================================= # Air Patch ========================== diff --git a/examples/2D_lid_driven_cavity/case.py b/examples/2D_lid_driven_cavity/case.py index 5e435f037..f206396be 100644 --- a/examples/2D_lid_driven_cavity/case.py +++ b/examples/2D_lid_driven_cavity/case.py @@ -44,6 +44,7 @@ 'bc_y%beg' : -16, 'bc_y%end' : -16, 'bc_y%ve1' : 0.5, + 'viscous' : 'T', # ========================================================== # Formatted Database Files Structure Parameters ============ @@ -79,4 +80,4 @@ 'fluid_pp(2)%Re(1)' : 1e4, # ========================================================== })) -# ============================================================================== \ No newline at end of file +# ============================================================================== diff --git a/examples/2D_mixing_artificial_Ma/case.py b/examples/2D_mixing_artificial_Ma/case.py index 4e28943bf..f1a6e8da0 100644 --- a/examples/2D_mixing_artificial_Ma/case.py +++ b/examples/2D_mixing_artificial_Ma/case.py @@ -96,6 +96,7 @@ 'bc_x%end' : -1, 'bc_y%beg' : -6, 'bc_y%end' : -6, + 'viscous' : 'T', # ========================================================================== # Formatted Database Files Structure Parameters ============================ diff --git a/examples/2D_rayleigh_taylor/case.py b/examples/2D_rayleigh_taylor/case.py index 23122d2c7..3f34263d3 100644 --- a/examples/2D_rayleigh_taylor/case.py +++ b/examples/2D_rayleigh_taylor/case.py @@ -69,6 +69,7 @@ 'bc_y%end' : -16, 'num_patches' : 1, 'num_fluids' : 2, + 'viscous' : 'T', # ======================================= # Database Structure Parameters ========= diff --git a/examples/2D_shearlayer/case.py b/examples/2D_shearlayer/case.py index 0474322d1..0c52cba96 100644 --- a/examples/2D_shearlayer/case.py +++ b/examples/2D_shearlayer/case.py @@ -45,6 +45,7 @@ 'bc_x%end' :-1, 'bc_y%beg' :-5, 'bc_y%end' :-5, + 'viscous' : 'T', # ========================================================================== # Formatted Database Files Structure Parameters ============================ diff --git a/examples/2D_viscous/case.py b/examples/2D_viscous/case.py index 99aabae0f..cb7058115 100644 --- a/examples/2D_viscous/case.py +++ b/examples/2D_viscous/case.py @@ -57,6 +57,7 @@ 'bc_x%end' : -1, 'bc_y%beg' : -6, 'bc_y%end' : -6, + 'viscous' : 'T', # ========================================================================== # Formatted Database Files Structure Parameters ============================ diff --git a/examples/3D_TaylorGreenVortex/case.py b/examples/3D_TaylorGreenVortex/case.py index 58c48b03c..31e116eb5 100644 --- a/examples/3D_TaylorGreenVortex/case.py +++ b/examples/3D_TaylorGreenVortex/case.py @@ -68,6 +68,7 @@ 'bc_y%end' : -1, 'bc_z%beg' : -1, 'bc_z%end' : -1, + 'viscous' : 'T', # ========================================================== # Formatted Database Files Structure Parameters ============ diff --git a/examples/3D_ibm_bowshock/case.py b/examples/3D_ibm_bowshock/case.py index a9f21a233..43e3eaa5d 100644 --- a/examples/3D_ibm_bowshock/case.py +++ b/examples/3D_ibm_bowshock/case.py @@ -70,6 +70,7 @@ # Set IB to True and add 1 patch 'ib' : 'T', 'num_ibs' : 1, + 'viscous' : 'T', # ========================================================================== # Formatted Database Files Structure Parameters ============================ @@ -113,4 +114,4 @@ 'fluid_pp(1)%pi_inf' : 0, 'fluid_pp(1)%Re(1)' : 7535533.2, # ========================================================================== -})) \ No newline at end of file +})) diff --git a/examples/3D_rayleigh_taylor/case.py b/examples/3D_rayleigh_taylor/case.py index a8b6b10d3..462eacc46 100644 --- a/examples/3D_rayleigh_taylor/case.py +++ b/examples/3D_rayleigh_taylor/case.py @@ -76,6 +76,7 @@ 'bc_z%end' : -3, 'num_patches' : 1, 'num_fluids' : 2, + 'viscous' : 'T', # ======================================= # Database Structure Parameters ========= diff --git a/examples/3D_recovering_sphere/case.py b/examples/3D_recovering_sphere/case.py index 3605bd3ff..68a2cad18 100644 --- a/examples/3D_recovering_sphere/case.py +++ b/examples/3D_recovering_sphere/case.py @@ -75,6 +75,7 @@ 'num_patches' : 2, 'num_fluids' : 2, 'weno_avg' : 'T', + 'surface_tension' : 'T', # ======================================= # Database Structure Parameters ========= @@ -91,13 +92,11 @@ # Fluid Parameters (Water) ============== 'fluid_pp(1)%gamma' : 1.E+00/(2.1E+00-1.E+00), 'fluid_pp(1)%pi_inf' : 2.1E+00*1.E+06/(2.1E+00-1.E+00), - #'fluid_pp(1)%Re(1)' : 1.e3, # ======================================= # Fluid Parameters (Gas) ================ 'fluid_pp(2)%gamma' : 1.E+00/(1.4E+00-1.E+00), 'fluid_pp(2)%pi_inf' : 0.E+00, - #'fluid_pp(2)%Re(1)' : 1.81e5, # ======================================= # Air Patch ========================== diff --git a/examples/3D_turb_mixing/case.py b/examples/3D_turb_mixing/case.py index cc9dd30a2..377b88fde 100644 --- a/examples/3D_turb_mixing/case.py +++ b/examples/3D_turb_mixing/case.py @@ -87,6 +87,7 @@ 'bc_y%end' : -6, 'bc_z%beg' : -1, 'bc_z%end' : -1, + 'viscous' : 'T', # ========================================================================== # Formatted Database Files Structure Parameters ============================ diff --git a/src/common/m_checker_common.fpp b/src/common/m_checker_common.fpp index c8a6fa33e..71ff6a185 100644 --- a/src/common/m_checker_common.fpp +++ b/src/common/m_checker_common.fpp @@ -290,11 +290,18 @@ contains !> Checks constraints on the surface tension parameters. !! Called by s_check_inputs_common for all three stages subroutine s_check_inputs_surface_tension - @:PROHIBIT(.not. f_is_default(sigma) .and. sigma < 0d0, & + @:PROHIBIT(surface_tension .and. sigma < 0d0, & "sigma must be greater than or equal to zero") - @:PROHIBIT(.not. f_is_default(sigma) .and. model_eqns /= 3, & + @:PROHIBIT(surface_tension .and. sigma == dflt_real, & + "sigma must be set if surface_tension is enabled") + + @:PROHIBIT(.not. f_is_default(sigma) .and. .not. surface_tension, & + "sigma is set but surface_tension is not enabled") + + @:PROHIBIT(surface_tension .and. model_eqns /= 3, & "The surface tension model requires model_eqns=3") + end subroutine s_check_inputs_surface_tension !> Checks constraints on the inputs for moving boundaries. diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index b02180f00..42748e57b 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -245,7 +245,7 @@ contains MPI_DOUBLE_PRECISION, MPI_MAX, 0, & MPI_COMM_WORLD, ierr) - if (any(Re_size > 0)) then + if (viscous) then call MPI_REDUCE(vcfl_max_loc, vcfl_max_glb, 1, & MPI_DOUBLE_PRECISION, MPI_MAX, 0, & MPI_COMM_WORLD, ierr) @@ -258,7 +258,7 @@ contains icfl_max_glb = icfl_max_loc - if (any(Re_size > 0)) then + if (viscous) then vcfl_max_glb = vcfl_max_loc Rc_min_glb = Rc_min_loc end if diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index c82a15c97..09462ce2c 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -342,7 +342,7 @@ contains #ifdef MFC_SIMULATION ! Computing the shear and bulk Reynolds numbers from species analogs - if (any(Re_size > 0)) then + if (viscous) then if (num_fluids == 1) then ! need to consider case with num_fluids >= 2 do i = 1, 2 @@ -532,7 +532,7 @@ contains G_K = max(0d0, G_K) end if - if (any(Re_size > 0)) then + if (viscous) then do i = 1, 2 Re_K(i) = dflt_real @@ -598,7 +598,7 @@ contains qv_K = qvs(1) end if - if (any(Re_size > 0)) then + if (viscous) then if (num_fluids == 1) then ! need to consider case with num_fluids >= 2 do i = 1, 2 @@ -663,7 +663,7 @@ contains #ifdef MFC_SIMULATION - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(Res(1:2, 1:maxval(Re_size))) do i = 1, 2 do j = 1, Re_size(i) @@ -1047,7 +1047,7 @@ contains qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l) end do - if (.not. f_is_default(sigma)) then + if (surface_tension) then qK_prim_vf(c_idx)%sf(j, k, l) = qK_cons_vf(c_idx)%sf(j, k, l) end if @@ -1221,7 +1221,7 @@ contains end do end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then q_cons_vf(c_idx)%sf(j, k, l) = q_prim_vf(c_idx)%sf(j, k, l) end if diff --git a/src/post_process/m_checker.fpp b/src/post_process/m_checker.fpp index 3ebdf3f87..2b2bf88e6 100644 --- a/src/post_process/m_checker.fpp +++ b/src/post_process/m_checker.fpp @@ -115,7 +115,7 @@ contains !> Checks constraints on surface tension parameters (cf_wrt and sigma) subroutine s_check_inputs_surface_tension - @:PROHIBIT(cf_wrt .and. f_is_default(sigma), & + @:PROHIBIT(cf_wrt .and. .not. surface_tension, & "cf_wrt can only be enabled if the surface coefficient is set") end subroutine s_check_inputs_surface_tension diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 260b27d3b..e1bf77379 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -275,6 +275,7 @@ module m_global_parameters !> @name surface tension coefficient !> @{ real(kind(0d0)) :: sigma + logical :: surface_tension !> #} !> @name Index variables used for m_variables_conversion @@ -397,6 +398,7 @@ contains poly_sigma = dflt_real sigR = dflt_real sigma = dflt_real + surface_tension = .false. adv_n = .false. end subroutine s_assign_default_values_to_user_inputs @@ -534,7 +536,7 @@ contains sys_size = stress_idx%end end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then c_idx = sys_size + 1 sys_size = c_idx end if @@ -559,7 +561,7 @@ contains sys_size = internalEnergies_idx%end alf_idx = 1 ! dummy, cannot actually have a void fraction - if (.not. f_is_default(sigma)) then + if (surface_tension) then c_idx = sys_size + 1 sys_size = c_idx end if diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90 index 90a2da971..7156f16db 100644 --- a/src/post_process/m_start_up.f90 +++ b/src/post_process/m_start_up.f90 @@ -82,7 +82,7 @@ subroutine s_read_input_file polydisperse, poly_sigma, file_per_process, relax, & relax_model, cf_wrt, sigma, adv_n, ib, & cfl_adap_dt, cfl_const_dt, t_save, t_stop, n_start, & - cfl_target + cfl_target, surface_tension ! Inquiring the status of the post_process.inp file file_loc = 'post_process.inp' diff --git a/src/pre_process/m_assign_variables.fpp b/src/pre_process/m_assign_variables.fpp index 8459f36a5..c24e0ee40 100644 --- a/src/pre_process/m_assign_variables.fpp +++ b/src/pre_process/m_assign_variables.fpp @@ -663,7 +663,7 @@ contains end do end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then q_prim_vf(c_idx)%sf(j, k, l) = eta*patch_icpp(patch_id)%cf_val + & (1d0 - eta)*patch_icpp(smooth_patch_id)%cf_val end if diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index c35698ebb..3d584b915 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -228,6 +228,7 @@ module m_global_parameters !> @name Surface Tension Modeling !> @{ real(kind(0d0)) :: sigma + logical :: surface_tension !> @} !> @name Index variables used for m_variables_conversion @@ -400,6 +401,7 @@ contains Re_inv = dflt_real Web = dflt_real poly_sigma = dflt_real + surface_tension = .false. adv_n = .false. @@ -615,7 +617,7 @@ contains sys_size = stress_idx%end end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then c_idx = sys_size + 1 sys_size = c_idx end if @@ -639,7 +641,7 @@ contains internalEnergies_idx%end = adv_idx%end + num_fluids sys_size = internalEnergies_idx%end - if (.not. f_is_default(sigma)) then + if (surface_tension) then c_idx = sys_size + 1 sys_size = c_idx end if diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp index 24558b098..33895c7c4 100644 --- a/src/pre_process/m_mpi_proxy.fpp +++ b/src/pre_process/m_mpi_proxy.fpp @@ -56,7 +56,7 @@ contains & 'perturb_flow', 'perturb_sph', 'mixlayer_vel_profile', & & 'mixlayer_perturb', 'bubbles', 'polytropic', 'polydisperse', & & 'qbmm', 'file_per_process', 'adv_n', 'ib' , 'cfl_adap_dt', & - & 'cfl_const_dt', 'cfl_dt'] + & 'cfl_const_dt', 'cfl_dt', 'surface_tension'] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor call MPI_BCAST(fluid_rho(1), num_fluids_max, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index 7712db0a0..641752644 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -141,7 +141,8 @@ contains sigR, sigV, dist_type, rhoRV, R0_type, & file_per_process, relax, relax_model, & palpha_eps, ptgalpha_eps, ib, num_ibs, patch_ib, & - sigma, adv_n, cfl_adap_dt, cfl_const_dt, n_start, n_start_old + sigma, adv_n, cfl_adap_dt, cfl_const_dt, n_start, & + n_start_old, surface_tension ! Inquiring the status of the pre_process.inp file file_loc = 'pre_process.inp' diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp index adaea9ab6..b45ac8f82 100644 --- a/src/simulation/m_checker.fpp +++ b/src/simulation/m_checker.fpp @@ -275,7 +275,14 @@ contains @:PROHIBIT(weno_order == 1 .and. (.not. weno_avg) .and. (.not. f_is_default(fluid_pp(i)%Re(j))), & "weno_order = 1 without weno_avg does not support fluid_pp("//trim(iStr)//")%"// "Re("//trim(jStr)//")") end do + @:PROHIBIT(.not. f_is_default(fluid_pp(i)%Re(1)) .and. .not. viscous, & + "Re(1) is specified, but viscous is not set to true") + @:PROHIBIT(.not. f_is_default(fluid_pp(i)%Re(2)) .and. .not. viscous, & + "Re(2) is specified, but viscous is not set to true") + @:PROHIBIT(f_is_default(fluid_pp(i)%Re(1)) .and. viscous, & + "Re(1) is not specified, but viscous is set to true") end do + end subroutine s_check_inputs_stiffened_eos_viscosity !> Checks constraints on body forces parameters (bf_x[y,z], etc.) diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 5e9386708..dcba87d60 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -148,7 +148,7 @@ contains ! Generating table header for the stability criteria to be outputted if (cfl_dt) then - if (any(Re_size > 0)) then + if (viscous) then write (1, '(A)') '==== Time-steps ====== dt ===== Time ======= ICFL '// & 'Max ==== VCFL Max ====== Rc Min =======' else @@ -156,7 +156,7 @@ contains '============== ICFL Max =============' end if else - if (any(Re_size > 0)) then + if (viscous) then write (1, '(A)') '==== Time-steps ====== Time ======= ICFL '// & 'Max ==== VCFL Max ====== Rc Min =======' else @@ -262,7 +262,7 @@ contains ! Compute mixture sound speed call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0d0, c) - if (any(Re_size > 0)) then + if (viscous) then call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf, vcfl_sf, Rc_sf) else call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf) @@ -278,13 +278,13 @@ contains #ifdef CRAY_ACC_WAR !$acc update host(icfl_sf) - if (any(Re_size > 0)) then + if (viscous) then !$acc update host(vcfl_sf, Rc_sf) end if icfl_max_loc = maxval(icfl_sf) - if (any(Re_size > 0)) then + if (viscous) then vcfl_max_loc = maxval(vcfl_sf) Rc_min_loc = minval(Rc_sf) end if @@ -293,7 +293,7 @@ contains icfl_max_loc = maxval(icfl_sf) !$acc end kernels - if (any(Re_size > 0)) then + if (viscous) then !$acc kernels vcfl_max_loc = maxval(vcfl_sf) Rc_min_loc = minval(Rc_sf) @@ -313,21 +313,21 @@ contains Rc_min_glb) else icfl_max_glb = icfl_max_loc - if (any(Re_size > 0)) vcfl_max_glb = vcfl_max_loc - if (any(Re_size > 0)) Rc_min_glb = Rc_min_loc + if (viscous) vcfl_max_glb = vcfl_max_loc + if (viscous) Rc_min_glb = Rc_min_loc end if ! Determining the stability criteria extrema over all the time-steps if (icfl_max_glb > icfl_max) icfl_max = icfl_max_glb - if (any(Re_size > 0)) then + if (viscous) then if (vcfl_max_glb > vcfl_max) vcfl_max = vcfl_max_glb if (Rc_min_glb < Rc_min) Rc_min = Rc_min_glb end if ! Outputting global stability criteria extrema at current time-step if (proc_rank == 0) then - if (any(Re_size > 0)) then + if (viscous) then write (1, '(6X,I8,F10.6,6X,6X,F10.6,6X,F9.6,6X,F9.6,6X,F10.6)') & t_step, dt, t_step*dt, icfl_max_glb, & vcfl_max_glb, & @@ -344,7 +344,7 @@ contains call s_mpi_abort('ICFL is greater than 1.0. Exiting ...') end if - if (any(Re_size > 0)) then + if (viscous) then if (vcfl_max_glb /= vcfl_max_glb) then call s_mpi_abort('VCFL is NaN. Exiting ...') elseif (vcfl_max_glb > 1d0) then @@ -1574,8 +1574,8 @@ contains write (3, '(A)') '' write (3, '(A,F9.6)') 'ICFL Max: ', icfl_max - if (any(Re_size > 0)) write (3, '(A,F9.6)') 'VCFL Max: ', vcfl_max - if (any(Re_size > 0)) write (3, '(A,F10.6)') 'Rc Min: ', Rc_min + if (viscous) write (3, '(A,F9.6)') 'VCFL Max: ', vcfl_max + if (viscous) write (3, '(A,F10.6)') 'Rc Min: ', Rc_min call cpu_time(run_time) @@ -1611,7 +1611,7 @@ contains @:ALLOCATE_GLOBAL(icfl_sf(0:m, 0:n, 0:p)) icfl_max = 0d0 - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(vcfl_sf(0:m, 0:n, 0:p)) @:ALLOCATE_GLOBAL(Rc_sf (0:m, 0:n, 0:p)) @@ -1648,7 +1648,7 @@ contains ! Deallocating the ICFL, VCFL, CCFL, and Rc stability criteria @:DEALLOCATE_GLOBAL(icfl_sf) - if (any(Re_size > 0)) then + if (viscous) then @:DEALLOCATE_GLOBAL(vcfl_sf, Rc_sf) end if diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index dd5808ebb..63fc21c73 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -163,6 +163,9 @@ module m_global_parameters logical :: hypoelasticity !< hypoelasticity modeling logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling logical :: cu_tensor + logical :: viscous !< Viscous effects + logical :: shear_stress !< Shear stresses + logical :: bulk_stress !< Bulk stresses !$acc declare create(chemistry) @@ -183,7 +186,7 @@ module m_global_parameters !$acc declare create(num_dims, weno_polyn, weno_order, weno_num_stencils, num_fluids, wenojs, mapped_weno, wenoz, teno, wenoz_q) #:endif - !$acc declare create(mpp_lim, model_eqns, mixture_err, alt_soundspeed, avg_state, mp_weno, weno_eps, teno_CT, hypoelasticity, low_Mach) + !$acc declare create(mpp_lim, model_eqns, mixture_err, alt_soundspeed, avg_state, mp_weno, weno_eps, teno_CT, hypoelasticity, low_Mach, viscous, shear_stress, bulk_stress) logical :: relax !< activate phase change integer :: relax_model !< Relaxation model @@ -459,7 +462,8 @@ module m_global_parameters !> @name Surface tension parameters !> @{ real(kind(0d0)) :: sigma - !$acc declare create(sigma) + logical :: surface_tension + !$acc declare create(sigma, surface_tension) !> @} integer :: momxb, momxe @@ -562,6 +566,9 @@ contains weno_flat = .true. riemann_flat = .true. rdma_mpi = .false. + viscous = .false. + shear_stress = .false. + bulk_stress = .false. #:if not MFC_CASE_OPTIMIZATION mapped_weno = .false. @@ -650,6 +657,7 @@ contains ! Surface tension sigma = dflt_real + surface_tension = .false. ! Cuda aware MPI cu_tensor = .false. @@ -888,7 +896,7 @@ contains sys_size = stress_idx%end end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then c_idx = sys_size + 1 sys_size = c_idx end if @@ -906,7 +914,7 @@ contains internalEnergies_idx%end = adv_idx%end + num_fluids sys_size = internalEnergies_idx%end - if (.not. f_is_default(sigma)) then + if (surface_tension) then c_idx = sys_size + 1 sys_size = c_idx end if @@ -973,11 +981,14 @@ contains if (fluid_pp(i)%Re(2) > 0) Re_size(2) = Re_size(2) + 1 end do - !$acc update device(Re_size) + if (Re_size(1) > 0d0) shear_stress = .true. + if (Re_size(2) > 0d0) bulk_stress = .true. + + !$acc update device(Re_size, viscous, shear_stress, bulk_stress) ! Bookkeeping the indexes of any viscous fluids and any pairs of ! fluids whose interface will support effects of surface tension - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(Re_idx(1:2, 1:maxval(Re_size))) @@ -1049,7 +1060,7 @@ contains ! sufficient boundary conditions data as to iterate the solution in ! the physical computational domain from one time-step iteration to ! the next one - if (any(Re_size > 0)) then + if (viscous) then buff_size = 2*weno_polyn + 2 ! else if (hypoelasticity) then !TODO: check if necessary ! buff_size = 2*weno_polyn + 2 @@ -1193,7 +1204,7 @@ contains ! Deallocating the variables bookkeeping the indexes of any viscous ! fluids and any pairs of fluids whose interfaces supported effects ! of surface tension - if (any(Re_size > 0)) then + if (viscous) then @:DEALLOCATE_GLOBAL(Re_idx) end if diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index 7da0b04d6..bc0552b1d 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -143,7 +143,7 @@ contains v_size = sys_size end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then nVars = num_dims + 1 if (n > 0) then if (p > 0) then @@ -203,7 +203,8 @@ contains & 'polydisperse', 'qbmm', 'acoustic_source', 'probe_wrt', 'integral_wrt', & & 'prim_vars_wrt', 'weno_avg', 'file_per_process', 'relax', & & 'adv_n', 'adap_dt', 'ib', 'bodyForces', 'bf_x', 'bf_y', 'bf_z', & - & 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt' ] + & 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt', 'surface_tension', & + & 'viscous', 'shear_stress', 'bulk_stress' ] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor @@ -2343,7 +2344,7 @@ contains @:DEALLOCATE_GLOBAL(ib_buff_send, ib_buff_recv) end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then @:DEALLOCATE_GLOBAL(c_divs_buff_send, c_divs_buff_recv) end if diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index a415baf21..084649837 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -246,7 +246,7 @@ contains @:ALLOCATE(q_prim_qp%vf(l)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end)) end do - if (.not. f_is_default(sigma)) then + if (surface_tension) then ! This assumes that the color function advection equation is ! the last equation. If this changes then this logic will ! need updated @@ -274,14 +274,14 @@ contains !$acc enter data attach(q_prim_qp%vf(l)%sf) end do - if (.not. f_is_default(sigma)) then + if (surface_tension) then q_prim_qp%vf(c_idx)%sf => & q_cons_qp%vf(c_idx)%sf !$acc enter data copyin(q_prim_qp%vf(c_idx)%sf) !$acc enter data attach(q_prim_qp%vf(c_idx)%sf) end if - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(tau_Re_vf(1:sys_size)) do i = 1, num_dims @:ALLOCATE(tau_Re_vf(cont_idx%end + i)%sf(idwbuff(1)%beg:idwbuff(1)%end, & @@ -377,43 +377,40 @@ contains @:ALLOCATE_GLOBAL(dq_prim_dy_qp(1:1)) @:ALLOCATE_GLOBAL(dq_prim_dz_qp(1:1)) - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE(dq_prim_dx_qp(1)%vf(1:sys_size)) @:ALLOCATE(dq_prim_dy_qp(1)%vf(1:sys_size)) @:ALLOCATE(dq_prim_dz_qp(1)%vf(1:sys_size)) - if (any(Re_size > 0)) then + + do l = mom_idx%beg, mom_idx%end + @:ALLOCATE(dq_prim_dx_qp(1)%vf(l)%sf( & + & idwbuff(1)%beg:idwbuff(1)%end, & + & idwbuff(2)%beg:idwbuff(2)%end, & + & idwbuff(3)%beg:idwbuff(3)%end)) + end do + + @:ACC_SETUP_VFs(dq_prim_dx_qp(1)) + + if (n > 0) then do l = mom_idx%beg, mom_idx%end - @:ALLOCATE(dq_prim_dx_qp(1)%vf(l)%sf( & - & idwbuff(1)%beg:idwbuff(1)%end, & - & idwbuff(2)%beg:idwbuff(2)%end, & - & idwbuff(3)%beg:idwbuff(3)%end)) + @:ALLOCATE(dq_prim_dy_qp(1)%vf(l)%sf( & + & idwbuff(1)%beg:idwbuff(1)%end, & + & idwbuff(2)%beg:idwbuff(2)%end, & + & idwbuff(3)%beg:idwbuff(3)%end)) end do - @:ACC_SETUP_VFs(dq_prim_dx_qp(1)) + @:ACC_SETUP_VFs(dq_prim_dy_qp(1)) - if (n > 0) then + if (p > 0) then do l = mom_idx%beg, mom_idx%end - @:ALLOCATE(dq_prim_dy_qp(1)%vf(l)%sf( & + @:ALLOCATE(dq_prim_dz_qp(1)%vf(l)%sf( & & idwbuff(1)%beg:idwbuff(1)%end, & & idwbuff(2)%beg:idwbuff(2)%end, & & idwbuff(3)%beg:idwbuff(3)%end)) end do - - @:ACC_SETUP_VFs(dq_prim_dy_qp(1)) - - if (p > 0) then - - do l = mom_idx%beg, mom_idx%end - @:ALLOCATE(dq_prim_dz_qp(1)%vf(l)%sf( & - & idwbuff(1)%beg:idwbuff(1)%end, & - & idwbuff(2)%beg:idwbuff(2)%end, & - & idwbuff(3)%beg:idwbuff(3)%end)) - end do - @:ACC_SETUP_VFs(dq_prim_dz_qp(1)) - end if - + @:ACC_SETUP_VFs(dq_prim_dz_qp(1)) end if end if @@ -446,7 +443,7 @@ contains @:ALLOCATE_GLOBAL(dqR_prim_dy_n(1:num_dims)) @:ALLOCATE_GLOBAL(dqR_prim_dz_n(1:num_dims)) - if (any(Re_size > 0)) then + if (viscous) then do i = 1, num_dims @:ALLOCATE(dqL_prim_dx_n(i)%vf(1:sys_size)) @:ALLOCATE(dqL_prim_dy_n(i)%vf(1:sys_size)) @@ -455,45 +452,41 @@ contains @:ALLOCATE(dqR_prim_dy_n(i)%vf(1:sys_size)) @:ALLOCATE(dqR_prim_dz_n(i)%vf(1:sys_size)) - if (any(Re_size > 0)) then + do l = mom_idx%beg, mom_idx%end + @:ALLOCATE(dqL_prim_dx_n(i)%vf(l)%sf( & + & idwbuff(1)%beg:idwbuff(1)%end, & + & idwbuff(2)%beg:idwbuff(2)%end, & + & idwbuff(3)%beg:idwbuff(3)%end)) + @:ALLOCATE(dqR_prim_dx_n(i)%vf(l)%sf( & + & idwbuff(1)%beg:idwbuff(1)%end, & + & idwbuff(2)%beg:idwbuff(2)%end, & + & idwbuff(3)%beg:idwbuff(3)%end)) + end do + if (n > 0) then do l = mom_idx%beg, mom_idx%end - @:ALLOCATE(dqL_prim_dx_n(i)%vf(l)%sf( & + @:ALLOCATE(dqL_prim_dy_n(i)%vf(l)%sf( & & idwbuff(1)%beg:idwbuff(1)%end, & & idwbuff(2)%beg:idwbuff(2)%end, & & idwbuff(3)%beg:idwbuff(3)%end)) - @:ALLOCATE(dqR_prim_dx_n(i)%vf(l)%sf( & + @:ALLOCATE(dqR_prim_dy_n(i)%vf(l)%sf( & & idwbuff(1)%beg:idwbuff(1)%end, & & idwbuff(2)%beg:idwbuff(2)%end, & & idwbuff(3)%beg:idwbuff(3)%end)) end do + end if - if (n > 0) then - do l = mom_idx%beg, mom_idx%end - @:ALLOCATE(dqL_prim_dy_n(i)%vf(l)%sf( & - & idwbuff(1)%beg:idwbuff(1)%end, & - & idwbuff(2)%beg:idwbuff(2)%end, & - & idwbuff(3)%beg:idwbuff(3)%end)) - @:ALLOCATE(dqR_prim_dy_n(i)%vf(l)%sf( & - & idwbuff(1)%beg:idwbuff(1)%end, & - & idwbuff(2)%beg:idwbuff(2)%end, & - & idwbuff(3)%beg:idwbuff(3)%end)) - end do - end if - - if (p > 0) then - do l = mom_idx%beg, mom_idx%end - @:ALLOCATE(dqL_prim_dz_n(i)%vf(l)%sf( & - & idwbuff(1)%beg:idwbuff(1)%end, & - & idwbuff(2)%beg:idwbuff(2)%end, & - & idwbuff(3)%beg:idwbuff(3)%end)) - @:ALLOCATE(dqR_prim_dz_n(i)%vf(l)%sf( & - & idwbuff(1)%beg:idwbuff(1)%end, & - & idwbuff(2)%beg:idwbuff(2)%end, & - & idwbuff(3)%beg:idwbuff(3)%end)) - end do - end if - + if (p > 0) then + do l = mom_idx%beg, mom_idx%end + @:ALLOCATE(dqL_prim_dz_n(i)%vf(l)%sf( & + & idwbuff(1)%beg:idwbuff(1)%end, & + & idwbuff(2)%beg:idwbuff(2)%end, & + & idwbuff(3)%beg:idwbuff(3)%end)) + @:ALLOCATE(dqR_prim_dz_n(i)%vf(l)%sf( & + & idwbuff(1)%beg:idwbuff(1)%end, & + & idwbuff(2)%beg:idwbuff(2)%end, & + & idwbuff(3)%beg:idwbuff(3)%end)) + end do end if @:ACC_SETUP_VFs(dqL_prim_dx_n(i), dqL_prim_dy_n(i), dqL_prim_dz_n(i)) @@ -502,7 +495,7 @@ contains end if ! END: Allocation/Association of d K_prim_ds_n ================== - if (any(Re_size > 0)) then + if (viscous) then if (weno_Re_flux) then @:ALLOCATE_GLOBAL(dqL_rsx_vf(idwbuff(1)%beg:idwbuff(1)%end, & idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, mom_idx%beg:mom_idx%end)) @@ -568,7 +561,7 @@ contains & idwbuff(3)%beg:idwbuff(3)%end)) end do - if (any(Re_size > 0) .or. (.not. f_is_default(sigma))) then + if (viscous .or. surface_tension) then do l = mom_idx%beg, E_idx @:ALLOCATE(flux_src_n(i)%vf(l)%sf( & & idwbuff(1)%beg:idwbuff(1)%end, & @@ -643,11 +636,11 @@ contains end do !$acc update device(gamma_min, pres_inf) - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(Res(1:2, 1:maxval(Re_size))) end if - if (any(Re_size > 0)) then + if (viscous) then do i = 1, 2 do j = 1, Re_size(i) Res(i, j) = fluid_pp(Re_idx(i, j))%Re(i) @@ -793,19 +786,19 @@ contains if (qbmm) call s_mom_inv(q_cons_qp%vf, q_prim_qp%vf, mom_sp, mom_3d, pb, rhs_pb, mv, rhs_mv, idwbuff(1), idwbuff(2), idwbuff(3), nbub) call nvtxStartRange("Viscous") - if (any(Re_size > 0)) call s_get_viscous(qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & - dqL_prim_dx_n, dqL_prim_dy_n, dqL_prim_dz_n, & - qL_prim, & - qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & - dqR_prim_dx_n, dqR_prim_dy_n, dqR_prim_dz_n, & - qR_prim, & - q_prim_qp, & - dq_prim_dx_qp, dq_prim_dy_qp, dq_prim_dz_qp, & - idwbuff(1), idwbuff(2), idwbuff(3)) + if (viscous) call s_get_viscous(qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & + dqL_prim_dx_n, dqL_prim_dy_n, dqL_prim_dz_n, & + qL_prim, & + qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & + dqR_prim_dx_n, dqR_prim_dy_n, dqR_prim_dz_n, & + qR_prim, & + q_prim_qp, & + dq_prim_dx_qp, dq_prim_dy_qp, dq_prim_dz_qp, & + idwbuff(1), idwbuff(2), idwbuff(3)) call nvtxEndRange call nvtxStartRange("Surface_Tension") - if (.not. f_is_default(sigma)) call s_get_capilary(q_prim_qp%vf) + if (surface_tension) call s_get_capilary(q_prim_qp%vf) call nvtxEndRange ! Dimensional Splitting Loop ======================================= @@ -815,7 +808,7 @@ contains call nvtxStartRange("RHS-WENO") - if (f_is_default(sigma)) then + if (.not. surface_tension) then ! Reconstruct densitiess iv%beg = 1; iv%end = sys_size call s_reconstruct_cell_boundary_values( & @@ -927,7 +920,7 @@ contains ! RHS additions for viscosity call nvtxStartRange("RHS_add_phys") - if (any(Re_size > 0d0) .or. (.not. f_is_default(sigma))) then + if (viscous .or. surface_tension) then call s_compute_additional_physics_rhs(id, & q_prim_qp%vf, & rhs_vf, & @@ -1604,7 +1597,7 @@ contains if (idir == 1) then ! x-direction - if (.not. f_is_default(sigma)) then + if (surface_tension) then !$acc parallel loop collapse(3) gang vector default(present) do l = 0, p do k = 0, n @@ -1636,7 +1629,7 @@ contains elseif (idir == 2) then ! y-direction - if (.not. f_is_default(sigma)) then + if (surface_tension) then !$acc parallel loop collapse(3) gang vector default(present) do l = 0, p do k = 0, n @@ -1652,7 +1645,7 @@ contains end if if (cyl_coord .and. ((bc_y%beg == -2) .or. (bc_y%beg == -14))) then - if (any(Re_size > 0)) then + if (viscous) then if (p > 0) then call s_compute_viscous_stress_tensor(q_prim_vf, & dq_prim_dx_vf(mom_idx%beg:mom_idx%end), & @@ -1736,7 +1729,7 @@ contains end do end do - if (any(Re_size > 0)) then + if (viscous) then !$acc parallel loop collapse(2) gang vector default(present) do l = 0, p do j = 0, m @@ -1771,7 +1764,7 @@ contains elseif (idir == 3) then ! z-direction - if (.not. f_is_default(sigma)) then + if (surface_tension) then !$acc parallel loop collapse(3) gang vector default(present) do l = 0, p do k = 0, n @@ -2035,7 +2028,7 @@ contains pi_inf = pi_inf + alpha(i)*pi_infs(i) end do - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re(i) = dflt_real @@ -2248,7 +2241,7 @@ contains @:DEALLOCATE_GLOBAL(qL_rsz_vf, qR_rsz_vf) end if - if (any(Re_size > 0) .and. weno_Re_flux) then + if (viscous .and. weno_Re_flux) then @:DEALLOCATE_GLOBAL(dqL_rsx_vf, dqR_rsx_vf) if (n > 0) then @@ -2265,7 +2258,7 @@ contains deallocate (alf_sum%sf) end if - if (any(Re_size > 0)) then + if (viscous) then do l = mom_idx%beg, mom_idx%end @:DEALLOCATE(dq_prim_dx_qp(1)%vf(l)%sf) end do @@ -2289,29 +2282,26 @@ contains @:DEALLOCATE(dq_prim_dz_qp(1)%vf) end if - if (any(Re_size > 0)) then + if (viscous) then do i = num_dims, 1, -1 - if (any(Re_size > 0)) then + do l = mom_idx%beg, mom_idx%end + @:DEALLOCATE(dqL_prim_dx_n(i)%vf(l)%sf) + @:DEALLOCATE(dqR_prim_dx_n(i)%vf(l)%sf) + end do + + if (n > 0) then do l = mom_idx%beg, mom_idx%end - @:DEALLOCATE(dqL_prim_dx_n(i)%vf(l)%sf) - @:DEALLOCATE(dqR_prim_dx_n(i)%vf(l)%sf) + @:DEALLOCATE(dqL_prim_dy_n(i)%vf(l)%sf) + @:DEALLOCATE(dqR_prim_dy_n(i)%vf(l)%sf) end do + end if - if (n > 0) then - do l = mom_idx%beg, mom_idx%end - @:DEALLOCATE(dqL_prim_dy_n(i)%vf(l)%sf) - @:DEALLOCATE(dqR_prim_dy_n(i)%vf(l)%sf) - end do - end if - - if (p > 0) then - do l = mom_idx%beg, mom_idx%end - @:DEALLOCATE(dqL_prim_dz_n(i)%vf(l)%sf) - @:DEALLOCATE(dqR_prim_dz_n(i)%vf(l)%sf) - end do - end if - + if (p > 0) then + do l = mom_idx%beg, mom_idx%end + @:DEALLOCATE(dqL_prim_dz_n(i)%vf(l)%sf) + @:DEALLOCATE(dqR_prim_dz_n(i)%vf(l)%sf) + end do end if @:DEALLOCATE(dqL_prim_dx_n(i)%vf) @@ -2339,7 +2329,7 @@ contains @:DEALLOCATE(flux_gsrc_n(i)%vf(l)%sf) end do - if (any(Re_size > 0)) then + if (viscous) then do l = mom_idx%beg, E_idx @:DEALLOCATE(flux_src_n(i)%vf(l)%sf) end do @@ -2363,7 +2353,7 @@ contains @:DEALLOCATE_GLOBAL(flux_n, flux_src_n, flux_gsrc_n) - if (any(Re_size > 0) .and. cyl_coord) then + if (viscous .and. cyl_coord) then do i = 1, num_dims @:DEALLOCATE(tau_re_vf(cont_idx%end + i)%sf) end do diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 9930e3ee6..99da073bb 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -456,7 +456,7 @@ contains qv_R = qv_R + alpha_rho_R(i)*qvs(i) end do - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_L(i) = dflt_real @@ -590,7 +590,7 @@ contains call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, & vel_avg_rms, c_sum_Yi_Phi, c_avg) - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_avg_rs${XYZ}$_vf(j, k, l, i) = 2d0/(1d0/Re_L(i) + 1d0/Re_R(i)) @@ -827,7 +827,7 @@ contains #:endfor - if (any(Re_size > 0)) then + if (viscous) then if (weno_Re_flux) then call s_compute_viscous_source_flux( & @@ -1085,7 +1085,7 @@ contains alpha_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, advxb + i - 1) end do - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_L(i) = dflt_real @@ -1139,7 +1139,7 @@ contains call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, & vel_avg_rms, 0d0, c_avg) - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_avg_rs${XYZ}$_vf(j, k, l, i) = 2d0/(1d0/Re_L(i) + 1d0/Re_R(i)) @@ -1205,7 +1205,7 @@ contains end do flux_rs${XYZ}$_vf(j, k, l, E_idx) = (E_L + pres_L)*vel_L(dir_idx(1)) - if (.not. f_is_default(sigma)) then + if (surface_tension) then flux_rs${XYZ}$_vf(j, k, l, c_idx) = & qL_prim_rs${XYZ}$_vf(j, k, l, c_idx)*s_S end if @@ -1239,7 +1239,7 @@ contains end do flux_rs${XYZ}$_vf(j, k, l, E_idx) = (E_R + pres_R)*vel_R(dir_idx(1)) - if (.not. f_is_default(sigma)) then + if (surface_tension) then flux_rs${XYZ}$_vf(j, k, l, c_idx) = & qR_prim_rs${XYZ}$_vf(j + 1, k, l, c_idx)*s_S end if @@ -1280,7 +1280,7 @@ contains end do flux_rs${XYZ}$_vf(j, k, l, E_idx) = (E_Star + p_Star)*s_S - if (.not. f_is_default(sigma)) then + if (surface_tension) then flux_rs${XYZ}$_vf(j, k, l, c_idx) = & qL_prim_rs${XYZ}$_vf(j, k, l, c_idx)*s_S end if @@ -1323,7 +1323,7 @@ contains ! Compute the star velocities for the non-conservative terms end do - if (.not. f_is_default(sigma)) then + if (surface_tension) then flux_rs${XYZ}$_vf(j, k, l, c_idx) = & qR_prim_rs${XYZ}$_vf(j + 1, k, l, c_idx)*s_S end if @@ -1691,7 +1691,7 @@ contains qv_R = qvs(1) end if - if (any(Re_size > 0)) then + if (viscous) then if (num_fluids == 1) then ! Need to consider case with num_fluids >= 2 !$acc loop seq do i = 1, 2 @@ -1854,7 +1854,7 @@ contains call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, & vel_avg_rms, 0d0, c_avg) - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_avg_rs${XYZ}$_vf(j, k, l, i) = 2d0/(1d0/Re_L(i) + 1d0/Re_R(i)) @@ -2157,7 +2157,7 @@ contains qv_R = qv_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i)*qvs(i) end do - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_L(i) = dflt_real @@ -2263,7 +2263,7 @@ contains call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, & vel_avg_rms, c_sum_Yi_Phi, c_avg) - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_avg_rs${XYZ}$_vf(j, k, l, i) = 2d0/(1d0/Re_L(i) + 1d0/Re_R(i)) @@ -2469,7 +2469,7 @@ contains #:endfor ! Computing HLLC flux and source flux for Euler system of equations - if (any(Re_size > 0)) then + if (viscous) then if (weno_Re_flux) then call s_compute_viscous_source_flux( & qL_prim_vf(momxb:momxe), & @@ -2495,7 +2495,7 @@ contains end if end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then call s_compute_capilary_source_flux( & q_prim_vf, & vel_src_rsx_vf, & @@ -2528,11 +2528,11 @@ contains end do !$acc update device(Gs) - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(Res(1:2, 1:maxval(Re_size))) end if - if (any(Re_size > 0)) then + if (viscous) then do i = 1, 2 do j = 1, Re_size(i) Res(i, j) = fluid_pp(Re_idx(i, j))%Re(i) @@ -2578,7 +2578,7 @@ contains @:ALLOCATE_GLOBAL(mom_sp_rsx_vf(is1%beg:is1%end + 1, is2%beg:is2%end, is3%beg:is3%end, 1:4)) end if - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(Re_avg_rsx_vf(is1%beg:is1%end, & is2%beg:is2%end, & is3%beg:is3%end, 1:2)) @@ -2606,7 +2606,7 @@ contains @:ALLOCATE_GLOBAL(mom_sp_rsy_vf(is1%beg:is1%end + 1, is2%beg:is2%end, is3%beg:is3%end, 1:4)) end if - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(Re_avg_rsy_vf(is1%beg:is1%end, & is2%beg:is2%end, & is3%beg:is3%end, 1:2)) @@ -2634,7 +2634,7 @@ contains @:ALLOCATE_GLOBAL(mom_sp_rsz_vf(is1%beg:is1%end + 1, is2%beg:is2%end, is3%beg:is3%end, 1:4)) end if - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(Re_avg_rsz_vf(is1%beg:is1%end, & is2%beg:is2%end, & is3%beg:is3%end, 1:2)) @@ -2733,7 +2733,7 @@ contains end do end do - if (any(Re_size > 0)) then + if (viscous) then !$acc parallel loop collapse(3) gang vector default(present) do i = momxb, momxe do l = isz%beg, isz%end @@ -2788,7 +2788,7 @@ contains end do end do - if (any(Re_size > 0)) then + if (viscous) then !$acc parallel loop collapse(3) gang vector default(present) do i = momxb, momxe @@ -2847,7 +2847,7 @@ contains end do end do - if (any(Re_size > 0)) then + if (viscous) then !$acc parallel loop collapse(3) gang vector default(present) do i = momxb, momxe @@ -2897,7 +2897,7 @@ contains end do end do - if (any(Re_size > 0)) then + if (viscous) then !$acc parallel loop collapse(3) gang vector default(present) do i = momxb, momxe @@ -2950,7 +2950,7 @@ contains end do end do - if (any(Re_size > 0)) then + if (viscous) then !$acc parallel loop collapse(3) gang vector default(present) do i = momxb, momxe do k = isy%beg, isy%end @@ -2994,7 +2994,7 @@ contains end do end do - if (any(Re_size > 0)) then + if (viscous) then !$acc parallel loop collapse(3) gang vector default(present) do i = momxb, momxe do k = isy%beg, isy%end @@ -3069,7 +3069,7 @@ contains if (norm_dir == 1) then - if (any(Re_size > 0) .or. (.not. f_is_default(sigma))) then + if (viscous .or. (surface_tension)) then !$acc parallel loop collapse(4) gang vector default(present) do i = momxb, E_idx @@ -3102,7 +3102,7 @@ contains ! Reshaping Inputted Data in y-direction =========================== elseif (norm_dir == 2) then - if (any(Re_size > 0) .or. (.not. f_is_default(sigma))) then + if (viscous .or. (surface_tension)) then !$acc parallel loop collapse(4) gang vector default(present) do i = momxb, E_idx do l = is3%beg, is3%end @@ -3133,7 +3133,7 @@ contains ! Reshaping Inputted Data in z-direction =========================== else - if (any(Re_size > 0) .or. (.not. f_is_default(sigma))) then + if (viscous .or. (surface_tension)) then !$acc parallel loop collapse(4) gang vector default(present) do i = momxb, E_idx do j = is1%beg, is1%end @@ -3224,7 +3224,7 @@ contains ! Viscous Stresses in z-direction ================================== if (norm_dir == 1) then - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3250,7 +3250,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3278,7 +3278,7 @@ contains if (n == 0) return - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dy, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3320,7 +3320,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dy, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3352,7 +3352,7 @@ contains if (p == 0) return - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3393,7 +3393,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private( avg_vel, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3423,7 +3423,7 @@ contains ! Viscous Stresses in r-direction ================================== elseif (norm_dir == 2) then - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dy, tau_Re) do l = isz%beg, isz%end @@ -3473,7 +3473,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dy, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3508,7 +3508,7 @@ contains if (p == 0) return - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dy, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3553,7 +3553,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3583,7 +3583,7 @@ contains ! Viscous Stresses in theta-direction ================================== else - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dy, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3648,7 +3648,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dy, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3748,7 +3748,7 @@ contains ! Viscous Stresses in x-direction ================================== if (norm_dir == 1) then - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3774,7 +3774,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3802,7 +3802,7 @@ contains if (n == 0) return - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(dvel_avg_dx, dvel_avg_dy, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3843,7 +3843,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dy, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3871,7 +3871,7 @@ contains if (p == 0) return - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3911,7 +3911,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3941,7 +3941,7 @@ contains ! Viscous Stresses in y-direction ================================== elseif (norm_dir == 2) then - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, dvel_avg_dy, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3986,7 +3986,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, dvel_avg_dy, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -4017,7 +4017,7 @@ contains if (p == 0) return - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dy, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -4058,7 +4058,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -4088,7 +4088,7 @@ contains ! Viscous Stresses in z-direction ================================== else - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, dvel_avg_dy, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -4145,7 +4145,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, dvel_avg_dy, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -4371,7 +4371,7 @@ contains ! to convert mixture or species variables to the mixture variables ! s_convert_to_mixture_variables => null() - if (Re_size(1) > 0) then + if (viscous) then @:DEALLOCATE_GLOBAL(Re_avg_rsx_vf) end if @:DEALLOCATE_GLOBAL(vel_src_rsx_vf) @@ -4384,7 +4384,7 @@ contains if (n == 0) return - if (Re_size(1) > 0) then + if (viscous) then @:DEALLOCATE_GLOBAL(Re_avg_rsy_vf) end if @:DEALLOCATE_GLOBAL(vel_src_rsy_vf) @@ -4397,7 +4397,7 @@ contains if (p == 0) return - if (Re_size(1) > 0) then + if (viscous) then @:DEALLOCATE_GLOBAL(Re_avg_rsz_vf) end if @:DEALLOCATE_GLOBAL(vel_src_rsz_vf) diff --git a/src/simulation/m_sim_helpers.f90 b/src/simulation/m_sim_helpers.f90 index 9b1f1ce21..a4af91281 100644 --- a/src/simulation/m_sim_helpers.f90 +++ b/src/simulation/m_sim_helpers.f90 @@ -118,7 +118,7 @@ subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl_sf, vcfl dz(l)/(abs(vel(3)) + c)) end if - if (any(Re_size > 0)) then + if (viscous) then if (grid_geometry == 3) then vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) & @@ -145,7 +145,7 @@ subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl_sf, vcfl icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), & dy(k)/(abs(vel(2)) + c)) - if (any(Re_size > 0)) then + if (viscous) then vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k))**2d0 @@ -159,7 +159,7 @@ subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl_sf, vcfl !1D icfl_sf(j, k, l) = (dt/dx(j))*(abs(vel(1)) + c) - if (any(Re_size > 0)) then + if (viscous) then vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/dx(j)**2d0 @@ -213,7 +213,7 @@ subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) dz(l)/(abs(vel(3)) + c)) end if - if (any(Re_size > 0)) then + if (viscous) then if (grid_geometry == 3) then vcfl_dt = cfl_target*(min(dx(j), dy(k), fltr_dtheta)**2d0) & /minval(1/(rho*Re_l)) @@ -228,7 +228,7 @@ subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), & dy(k)/(abs(vel(2)) + c)) - if (any(Re_size > 0)) then + if (viscous) then vcfl_dt = cfl_target*(min(dx(j), dy(k))**2d0)/maxval((1/Re_l)/rho) end if @@ -236,7 +236,7 @@ subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) !1D icfl_dt = cfl_target*(dx(j)/(abs(vel(1)) + c)) - if (any(Re_size > 0)) then + if (viscous) then vcfl_dt = cfl_target*(dx(j)**2d0)/minval(1/(rho*Re_l)) end if diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index f40429b53..633f0e8ce 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -164,7 +164,8 @@ contains pi_fac, adv_n, adap_dt, bf_x, bf_y, bf_z, & k_x, k_y, k_z, w_x, w_y, w_z, p_x, p_y, p_z, & g_x, g_y, g_z, n_start, t_save, t_stop, & - cfl_adap_dt, cfl_const_dt, cfl_target + cfl_adap_dt, cfl_const_dt, cfl_target, & + viscous, surface_tension ! Checking that an input file has been provided by the user. If it ! has, then the input file is read in, otherwise, simulation exits. @@ -1325,13 +1326,13 @@ contains call s_initialize_acoustic_src() end if - if (any(Re_size > 0)) then + if (viscous) then call s_initialize_viscous_module() end if call s_initialize_rhs_module() - if (.not. f_is_default(sigma)) call s_initialize_surface_tension_module() + if (surface_tension) call s_initialize_surface_tension_module() #if defined(MFC_OpenACC) && defined(MFC_MEMORY_DUMP) call acc_present_dump() @@ -1470,7 +1471,7 @@ contains !$acc update device(R_n, R_v, phi_vn, phi_nv, Pe_c, Tw, pv, M_n, M_v, k_n, k_v, pb0, mass_n0, mass_v0, Pe_T, Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c, omegaN , mul0, ss, gamma_v, mu_v, gamma_m, gamma_n, mu_n, gam) !$acc update device(acoustic_source, num_source) - !$acc update device(sigma) + !$acc update device(sigma, surface_tension) !$acc update device(dx, dy, dz, x_cb, x_cc, y_cb, y_cc, z_cb, z_cc) @@ -1507,11 +1508,11 @@ contains call s_finalize_mpi_proxy_module() call s_finalize_global_parameters_module() if (relax) call s_finalize_relaxation_solver_module() - if (any(Re_size > 0)) then + if (viscous) then call s_finalize_viscous_module() end if - if (.not. f_is_default(sigma)) call s_finalize_surface_tension_module() + if (surface_tension) call s_finalize_surface_tension_module() if (bodyForces) call s_finalize_body_forces_module() ! Terminating MPI execution environment diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 4f78bb28b..e9b260370 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -189,7 +189,7 @@ contains end do end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then @:ALLOCATE(q_prim_vf(c_idx)%sf(idwbuff(1)%beg:idwbuff(1)%end, & idwbuff(2)%beg:idwbuff(2)%end, & idwbuff(3)%beg:idwbuff(3)%end)) diff --git a/src/simulation/m_viscous.fpp b/src/simulation/m_viscous.fpp index 139fccc7c..3e2c3bf70 100644 --- a/src/simulation/m_viscous.fpp +++ b/src/simulation/m_viscous.fpp @@ -95,7 +95,7 @@ contains end do end do end do - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(alpha_visc, alpha_rho_visc, Re_visc, tau_Re ) do l = is3_viscous%beg, is3_viscous%end do k = -1, 1 @@ -161,7 +161,7 @@ contains pi_inf_visc = pi_inf_visc + alpha_visc(i)*pi_infs(i) end do - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_visc(i) = dflt_real @@ -202,7 +202,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private(alpha_visc, alpha_rho_visc, Re_visc, tau_Re ) do l = is3_viscous%beg, is3_viscous%end do k = -1, 1 @@ -268,7 +268,7 @@ contains pi_inf_visc = pi_inf_visc + alpha_visc(i)*pi_infs(i) end do - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_visc(i) = dflt_real @@ -306,7 +306,7 @@ contains if (p == 0) return - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(alpha_visc, alpha_rho_visc, Re_visc, tau_Re ) do l = is3_viscous%beg, is3_viscous%end do k = -1, 1 @@ -372,7 +372,7 @@ contains pi_inf_visc = pi_inf_visc + alpha_visc(i)*pi_infs(i) end do - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_visc(i) = dflt_real @@ -414,7 +414,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private(alpha_visc, alpha_rho_visc, Re_visc, tau_Re ) do l = is3_viscous%beg, is3_viscous%end do k = -1, 1 @@ -480,7 +480,7 @@ contains pi_inf_visc = pi_inf_visc + alpha_visc(i)*pi_infs(i) end do - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_visc(i) = dflt_real @@ -1023,7 +1023,7 @@ contains is1_viscous, is2_viscous, is3_viscous) end if - if (any(Re_size > 0)) then + if (viscous) then if (weno_Re_flux) then if (norm_dir == 2) then !$acc parallel loop collapse(4) gang vector default(present) @@ -1124,7 +1124,7 @@ contains is1_viscous, is2_viscous, is3_viscous) end if - if (any(Re_size > 0)) then + if (viscous) then if (weno_Re_flux) then if (norm_dir == 2) then !$acc parallel loop collapse(4) gang vector default(present) diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index 7738ca65d..49be5732b 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -87,7 +87,8 @@ def analytic(self): 'num_ibs': ParamType.INT, 'cfl_dt': ParamType.LOG, 'n_start': ParamType.INT, - 'n_start_old': ParamType.INT + 'n_start_old': ParamType.INT, + 'surface_tension': ParamType.LOG, }) for ib_id in range(1, 10+1): @@ -225,6 +226,8 @@ def analytic(self): 't_save': ParamType.REAL, 'cfl_target': ParamType.REAL, 'low_Mach': ParamType.INT, + 'surface_tension': ParamType.LOG, + 'viscous': ParamType.LOG, }) for var in [ 'diffusion', 'reactions' ]: @@ -336,6 +339,7 @@ def analytic(self): 't_save': ParamType.REAL, 't_stop': ParamType.REAL, 'n_start': ParamType.INT, + 'surface_tension': ParamType.LOG, }) for cmp_id in range(1,3+1): diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py index cf7208a01..99c1795e1 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -75,7 +75,8 @@ def alter_bcs(dimInfo): cases.append(define_case_d(stack, f"bc={bc}", get_bc_mods(bc, dimInfo))) def alter_capillary(): - stack.push('', {'patch_icpp(1)%cf_val':1, 'patch_icpp(2)%cf_val':0, 'patch_icpp(3)%cf_val':1, 'sigma':1, 'model_eqns':3}) + stack.push('', {'patch_icpp(1)%cf_val':1, 'patch_icpp(2)%cf_val':0, 'patch_icpp(3)%cf_val':1, + 'sigma':1, 'model_eqns':3, 'surface_tension': 'T'}) cases.append(define_case_d(stack, [f"capillary=T","model_eqns=3"],{})) stack.pop() @@ -157,7 +158,8 @@ def alter_num_fluids(dimInfo): if num_fluids == 1: stack.push("Viscous", { - 'fluid_pp(1)%Re(1)' : 0.0001, 'dt' : 1e-11, 'patch_icpp(1)%vel(1)': 1.0}) + 'fluid_pp(1)%Re(1)' : 0.0001, 'dt' : 1e-11, 'patch_icpp(1)%vel(1)': 1.0, + 'viscous': 'T'}) alter_ib(dimInfo, six_eqn_model=True) @@ -175,7 +177,7 @@ def alter_num_fluids(dimInfo): stack.push("Viscous", { 'fluid_pp(1)%Re(1)' : 0.001, 'fluid_pp(1)%Re(2)' : 0.001, 'fluid_pp(2)%Re(1)' : 0.001, 'fluid_pp(2)%Re(2)' : 0.001, 'dt' : 1e-11, - 'patch_icpp(1)%vel(1)': 1.0}) + 'patch_icpp(1)%vel(1)': 1.0, 'viscous': 'T'}) alter_ib(dimInfo, six_eqn_model=True) @@ -206,7 +208,8 @@ def alter_2d(): stack.push("Viscous", { 'fluid_pp(1)%Re(1)' : 0.0001, 'fluid_pp(1)%Re(2)' : 0.0001, - 'fluid_pp(2)%Re(1)' : 0.0001, 'fluid_pp(2)%Re(2)' : 0.0001, 'dt' : 1e-11}) + 'fluid_pp(2)%Re(1)' : 0.0001, 'fluid_pp(2)%Re(2)' : 0.0001, 'dt' : 1e-11, + 'viscous': 'T'}) cases.append(define_case_d(stack, "", {'weno_Re_flux': 'F'})) cases.append(define_case_d(stack, "weno_Re_flux", {'weno_Re_flux': 'T'})) @@ -249,7 +252,8 @@ def alter_3d(): stack.push("Viscous", { 'fluid_pp(1)%Re(1)' : 0.0001, 'fluid_pp(1)%Re(2)' : 0.0001, - 'fluid_pp(2)%Re(1)' : 0.0001, 'fluid_pp(2)%Re(2)' : 0.0001, 'dt' : 1e-11 + 'fluid_pp(2)%Re(1)' : 0.0001, 'fluid_pp(2)%Re(2)' : 0.0001, 'dt' : 1e-11, + 'viscous': 'T' }) cases.append(define_case_d(stack, "", {'weno_Re_flux': 'F'})) @@ -533,7 +537,7 @@ def alter_instability_wave(dimInfo): 'mixlayer_vel_profile': 'T', 'mixlayer_domain': 1.475, 'mixlayer_vel_coef': 0.6, 'mixlayer_perturb': 'T', 'weno_Re_flux': 'T', 'weno_avg': 'T', 'mapped_weno': 'T', 'fluid_pp(1)%gamma': 0.16393442623, 'fluid_pp(1)%pi_inf': 22.312399959394575, - 'fluid_pp(1)%Re(1)': 1.6881644098979287, + 'fluid_pp(1)%Re(1)': 1.6881644098979287, 'viscous': 'T', 'patch_icpp(1)%x_centroid': 180.0, 'patch_icpp(1)%length_x': 360.0, 'patch_icpp(1)%y_centroid': 0.0, 'patch_icpp(1)%length_y': 360.0, 'patch_icpp(1)%vel(1)': 1.1966855884162177, 'patch_icpp(1)%vel(2)': 0.0, 'patch_icpp(1)%pres': 1.0, @@ -648,7 +652,7 @@ def alter_viscosity(dimInfo): # Viscosity & bubbles checks if len(dimInfo[0]) > 0: stack.push("Viscosity -> Bubbles", - {"fluid_pp(1)%Re(1)": 50, "bubbles": 'T'}) + {"fluid_pp(1)%Re(1)": 50, "bubbles": 'T', "viscous": 'T'}) stack.push('', { 'nb' : 1, 'fluid_pp(1)%gamma' : 0.16, 'fluid_pp(1)%pi_inf': 3515.0, From 39e410a4924664ce56d47d58cf9c320eaab20f5f Mon Sep 17 00:00:00 2001 From: Max Hawkins <37495064+max-Hawkins@users.noreply.github.com> Date: Thu, 7 Nov 2024 14:17:50 -0500 Subject: [PATCH 2/9] Specify CUDA version for login node builds (#697) Co-authored-by: Spencer Bryngelson Co-authored-by: Spencer Bryngelson --- CMakeLists.txt | 14 +++++++------- toolchain/bootstrap/modules.sh | 4 ++-- toolchain/modules | 2 ++ 3 files changed, 11 insertions(+), 9 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index c0bb7616b..22268e250 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -127,7 +127,7 @@ if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") if (CMAKE_BUILD_TYPE STREQUAL "Debug") add_compile_options( - -Wall + -Wall -fcheck=all,no-array-temps -fbacktrace -fimplicit-none @@ -142,7 +142,7 @@ if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") add_compile_options( $<$:-fallow-invalid-boz> $<$:-fallow-argument-mismatch> - $<$:-fcheck=bounds> + $<$:-fcheck=bounds> ) endif() elseif (CMAKE_Fortran_COMPILER_ID STREQUAL "Cray") @@ -155,7 +155,7 @@ elseif (CMAKE_Fortran_COMPILER_ID STREQUAL "Cray") ) add_link_options("SHELL:-hkeepfiles") - + if (CMAKE_BUILD_TYPE STREQUAL "Debug") add_compile_options( "SHELL:-h acc_model=auto_async_none" @@ -197,7 +197,7 @@ elseif ((CMAKE_Fortran_COMPILER_ID STREQUAL "NVHPC") OR (CMAKE_Fortran_COMPILER_ if (DEFINED ENV{MFC_CUDA_CC}) string(REGEX MATCHALL "[0-9]+" MFC_CUDA_CC $ENV{MFC_CUDA_CC}) - message(STATUS "Found $MFC_CUDA_CC specified. GPU code will be generated for ${MFC_CUDA_CC}.") + message(STATUS "Found $MFC_CUDA_CC specified. GPU code will be generated for compute capability(ies) ${MFC_CUDA_CC}.") endif() endif() @@ -231,7 +231,7 @@ if (CMAKE_BUILD_TYPE STREQUAL "Release") else() message(STATUS "IPO / LTO is NOT available") endif() - endif() + endif() endif() if (CMAKE_BUILD_TYPE STREQUAL "Debug") @@ -295,7 +295,7 @@ macro(HANDLE_SOURCES target useCommon) string(TOUPPER ${target} ${target}_UPPER) - # Gather: + # Gather: # * src/[,(common)]/*.f90 # * (if any) /modules//*.f90 file(GLOB ${target}_F90s CONFIGURE_DEPENDS "${${target}_DIR}/*.f90" @@ -388,7 +388,7 @@ function(MFC_SETUP_TARGET) # A little hacky, but it *is* an edge-case for *one* compiler. if (NVHPC_USE_TWO_PASS_IPO) add_library(${ARGS_TARGET}_lib OBJECT ${ARGS_SOURCES}) - target_compile_options(${ARGS_TARGET}_lib PRIVATE + target_compile_options(${ARGS_TARGET}_lib PRIVATE $<$:-Mextract=lib:${ARGS_TARGET}_lib> $<$:-Minline> ) diff --git a/toolchain/bootstrap/modules.sh b/toolchain/bootstrap/modules.sh index 3c83d9fdc..db964c436 100644 --- a/toolchain/bootstrap/modules.sh +++ b/toolchain/bootstrap/modules.sh @@ -91,8 +91,8 @@ if ! module load $MODULES; then fi if [ $(echo "$VARIABLES" | grep = | wc -c) -gt 0 ]; then - log " $ export $VARIABLES" - export $VARIABLES > /dev/null + log " $ export $(eval "echo $VARIABLES")" + export $(eval "echo $VARIABLES") > /dev/null fi # Don't check for Cray paths on Carpenter, otherwise do check if they exist diff --git a/toolchain/modules b/toolchain/modules index a2c558df4..2d0ebe925 100644 --- a/toolchain/modules +++ b/toolchain/modules @@ -45,6 +45,7 @@ p GT Phoenix p-all python/3.10.10 p-cpu gcc/12.3.0 openmpi/4.1.5 p-gpu nvhpc/24.5 hpcx/2.19-cuda cuda/12.1.1 +p-gpu MFC_CUDA_CC=70,80,89,90 NVHPC_CUDA_HOME=$CUDA_HOME CC=nvc CXX=nvc++ FC=nvfortran f OLCF Frontier f-all cce/18.0.0 cpe/24.07 rocm/6.1.3 cray-mpich/8.1.28 @@ -56,6 +57,7 @@ d-all python/3.11.6 d-cpu gcc/11.4.0 openmpi d-gpu nvhpc/24.1 cuda/12.3.0 openmpi/4.1.5+cuda cmake d-gpu CC=nvc CXX=nvc++ FC=nvfortran +d-gpu MFC_CUDA_CC=80,86 c DoD Carpenter c-all python/3.12.1 From 4396130ad71cb446691458c2ae78e8e0456278c9 Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Thu, 7 Nov 2024 20:43:10 -0500 Subject: [PATCH 3/9] Fix 2D_shockdroplet example case (#705) Co-authored-by: Spencer Bryngelson --- examples/2D_shockdroplet/case.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/2D_shockdroplet/case.py b/examples/2D_shockdroplet/case.py index 20df1d032..8a225deee 100755 --- a/examples/2D_shockdroplet/case.py +++ b/examples/2D_shockdroplet/case.py @@ -86,7 +86,7 @@ 'patch_icpp(1)%x_centroid' : 8*D, 'patch_icpp(1)%y_centroid' : 6*D, 'patch_icpp(1)%length_x' : 24*D, - 'patch_icpp(1)%length_y' : 12*D, + 'patch_icpp(1)%length_y' : 14*D, 'patch_icpp(1)%vel(1)' : 0., 'patch_icpp(1)%vel(2)' : 0.E+00, 'patch_icpp(1)%pres' : 101325., @@ -102,7 +102,7 @@ 'patch_icpp(2)%x_centroid' : -2.5*D, 'patch_icpp(2)%y_centroid' : 6*D, 'patch_icpp(2)%length_x' : 3*D, - 'patch_icpp(2)%length_y' : 12*D, + 'patch_icpp(2)%length_y' : 14*D, 'patch_icpp(2)%vel(1)' : vel, 'patch_icpp(2)%vel(2)' : 0.E+00, 'patch_icpp(2)%pres' : ps, @@ -134,4 +134,4 @@ 'fluid_pp(2)%gamma' : 1.E+00/(gam_a-1.E+00), 'fluid_pp(2)%pi_inf' : 0.E+00, # ========================================================== -})) \ No newline at end of file +})) From acd2763fa0feb7515bc6b5591d6bf1b5d19bb4b1 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Thu, 7 Nov 2024 23:52:23 -0500 Subject: [PATCH 4/9] Try to fix codecov (#707) --- .github/workflows/{coverage.yml => codecov.yml} | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) rename .github/workflows/{coverage.yml => codecov.yml} (91%) diff --git a/.github/workflows/coverage.yml b/.github/workflows/codecov.yml similarity index 91% rename from .github/workflows/coverage.yml rename to .github/workflows/codecov.yml index e779919b4..685215646 100644 --- a/.github/workflows/coverage.yml +++ b/.github/workflows/codecov.yml @@ -31,8 +31,7 @@ jobs: run: | sudo apt update -y sudo apt install -y tar wget make cmake gcc g++ python3 python3-dev "openmpi-*" libopenmpi-dev - - + - name: Build run: /bin/bash mfc.sh build -j $(nproc) --gcov @@ -40,10 +39,10 @@ jobs: run: /bin/bash mfc.sh test -a -j $(nproc) - name: Upload coverage reports to Codecov - uses: codecov/codecov-action@v4 + uses: codecov/codecov-action@v4.6.0 with: - token: ${{ secrets.CODECOV_TOKEN }} fail_ci_if_error: false + name: mfc-coverage verbose: true env: CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} From 82b53840aafe843c19fbdc37f8638a9224c21bd6 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 8 Nov 2024 00:25:50 -0500 Subject: [PATCH 5/9] Fix codecov maybe (#708) --- .github/workflows/codecov.yml | 55 ++++-------------------------- .github/workflows/coverage.yml | 48 ++++++++++++++++++++++++++ src/post_process/m_data_output.fpp | 2 +- 3 files changed, 55 insertions(+), 50 deletions(-) create mode 100644 .github/workflows/coverage.yml diff --git a/.github/workflows/codecov.yml b/.github/workflows/codecov.yml index 685215646..61a42ae96 100644 --- a/.github/workflows/codecov.yml +++ b/.github/workflows/codecov.yml @@ -1,49 +1,6 @@ -name: Coverage Check - -on: [push, pull_request, workflow_dispatch] - -jobs: - file-changes: - name: Detect File Changes - runs-on: 'ubuntu-latest' - outputs: - checkall: ${{ steps.changes.outputs.checkall }} - steps: - - name: Clone - uses: actions/checkout@v4 - - - name: Detect Changes - uses: dorny/paths-filter@v3 - id: changes - with: - filters: ".github/file-filter.yml" - - run: - name: Coverage Test on CodeCov - if: needs.file-changes.outputs.checkall == 'true' - needs: file-changes - runs-on: "ubuntu-latest" - steps: - - name: Checkouts - uses: actions/checkout@v4 - - - name: Setup Ubuntu - run: | - sudo apt update -y - sudo apt install -y tar wget make cmake gcc g++ python3 python3-dev "openmpi-*" libopenmpi-dev - - - name: Build - run: /bin/bash mfc.sh build -j $(nproc) --gcov - - - name: Test - run: /bin/bash mfc.sh test -a -j $(nproc) - - - name: Upload coverage reports to Codecov - uses: codecov/codecov-action@v4.6.0 - with: - fail_ci_if_error: false - name: mfc-coverage - verbose: true - env: - CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} - +coverage: + status: + project: + default: + target: 1% + threshold: 1% diff --git a/.github/workflows/coverage.yml b/.github/workflows/coverage.yml new file mode 100644 index 000000000..7487d8e55 --- /dev/null +++ b/.github/workflows/coverage.yml @@ -0,0 +1,48 @@ +name: Coverage Check + +on: [push, pull_request, workflow_dispatch] + +jobs: + file-changes: + name: Detect File Changes + runs-on: 'ubuntu-latest' + outputs: + checkall: ${{ steps.changes.outputs.checkall }} + steps: + - name: Clone + uses: actions/checkout@v4 + + - name: Detect Changes + uses: dorny/paths-filter@v3 + id: changes + with: + filters: ".github/file-filter.yml" + + run: + name: Coverage Test on CodeCov + if: needs.file-changes.outputs.checkall == 'true' + needs: file-changes + runs-on: "ubuntu-latest" + steps: + - name: Checkouts + uses: actions/checkout@v4 + + - name: Setup Ubuntu + run: | + sudo apt update -y + sudo apt install -y tar wget make cmake gcc g++ python3 python3-dev "openmpi-*" libopenmpi-dev + + - name: Build + run: /bin/bash mfc.sh build -j $(nproc) --gcov + + - name: Test + run: /bin/bash mfc.sh test -a -j $(nproc) + + - name: Upload coverage reports to Codecov + uses: codecov/codecov-action@v4 + with: + fail_ci_if_error: false + verbose: true + env: + CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} + diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp index 7aa7e8739..b29e5ab01 100644 --- a/src/post_process/m_data_output.fpp +++ b/src/post_process/m_data_output.fpp @@ -113,7 +113,7 @@ contains logical :: dir_check ! Generic loop iterator - integer :: i + integer :: i, k ,l ! Allocating the generic storage for the flow variable(s) that are ! going to be written to the formatted database file(s). Note once From 73782f55b566bc10aea1cfecc6aaabbb7e080e46 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 8 Nov 2024 09:37:36 -0500 Subject: [PATCH 6/9] Fixest format CI failure (#709) --- src/post_process/m_data_output.fpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp index b29e5ab01..784dd7c00 100644 --- a/src/post_process/m_data_output.fpp +++ b/src/post_process/m_data_output.fpp @@ -112,8 +112,7 @@ contains ! Generic logical used to test the existence of a particular folder logical :: dir_check - ! Generic loop iterator - integer :: i, k ,l + integer :: i ! Allocating the generic storage for the flow variable(s) that are ! going to be written to the formatted database file(s). Note once From 8775a6151b01fde2d4c8ad3e6e1e84a526449da6 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 8 Nov 2024 09:40:55 -0500 Subject: [PATCH 7/9] Remove [[some!]] unused variables (#699) --- src/common/m_phase_change.fpp | 2 -- src/simulation/m_bubbles.fpp | 14 ++------- src/simulation/m_chemistry.fpp | 11 ------- src/simulation/m_data_output.fpp | 43 ++++------------------------ src/simulation/m_mpi_proxy.fpp | 3 -- src/simulation/m_qbmm.fpp | 13 ++++----- src/simulation/m_rhs.fpp | 30 +++---------------- src/simulation/m_riemann_solvers.fpp | 15 ++-------- src/simulation/m_sim_helpers.f90 | 2 +- src/simulation/m_start_up.fpp | 7 +---- src/simulation/m_time_steppers.fpp | 16 ++--------- src/simulation/m_viscous.fpp | 2 -- src/simulation/m_weno.fpp | 6 +--- 13 files changed, 26 insertions(+), 138 deletions(-) diff --git a/src/common/m_phase_change.fpp b/src/common/m_phase_change.fpp index 9ae303ff3..cac791f70 100644 --- a/src/common/m_phase_change.fpp +++ b/src/common/m_phase_change.fpp @@ -313,8 +313,6 @@ contains real(kind(0.0d0)), intent(in) :: rhoe real(kind(0.0d0)), intent(out) :: TS - integer, dimension(num_fluids) :: ig !< flags to toggle the inclusion of fluids for the pT-equilibrium - real(kind(0.0d0)), dimension(num_fluids) :: pk !< individual initial pressures real(kind(0.0d0)) :: gp, gpp, hp, pO, mCP, mQ !< variables for the Newton Solver integer :: i, ns !< generic loop iterators diff --git a/src/simulation/m_bubbles.fpp b/src/simulation/m_bubbles.fpp index 804b473b8..82ec9e4e4 100644 --- a/src/simulation/m_bubbles.fpp +++ b/src/simulation/m_bubbles.fpp @@ -57,7 +57,7 @@ contains subroutine s_initialize_bubbles_module - integer :: i, j, k, l, q + integer :: l @:ALLOCATE_GLOBAL(rs(1:nb)) @:ALLOCATE_GLOBAL(vs(1:nb)) @@ -119,7 +119,7 @@ contains integer, intent(in) :: idir type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf - integer :: i, j, k, l, q + integer :: j, k, l if (idir == 1) then @@ -183,18 +183,11 @@ contains real(kind(0d0)) :: rddot real(kind(0d0)) :: pb, mv, vflux, pbdot real(kind(0d0)) :: n_tait, B_tait - real(kind(0d0)), dimension(nb) :: Rtmp, Vtmp real(kind(0d0)) :: myR, myV, alf, myP, myRho, R2Vav, R3 real(kind(0d0)), dimension(num_fluids) :: myalpha, myalpha_rho - real(kind(0d0)) :: start, finish - real(kind(0d0)) :: nbub !< Bubble number density - - real(kind(0d0)), dimension(2) :: Re !< Reynolds number - integer :: i, j, k, l, q, ii !< Loop variables - integer :: ndirs !< Number of coordinate directions real(kind(0d0)) :: err1, err2, err3, err4, err5 !< Error estimates for adaptive time stepping real(kind(0d0)) :: t_new !< Updated time step size @@ -452,7 +445,7 @@ contains real(kind(0d0)), intent(IN) :: fntait, fBtait, f_bub_adv_src, f_divu real(kind(0d0)), intent(out) :: h - real(kind(0d0)) :: h0, h1, h_min !< Time step size + real(kind(0d0)) :: h0, h1 !< Time step size real(kind(0d0)) :: d0, d1, d2 !< norms real(kind(0d0)), dimension(2) :: myR_tmp, myV_tmp, myA_tmp !< Bubble radius, radial velocity, and radial acceleration @@ -912,7 +905,6 @@ contains real(kind(0.d0)) :: T_bar real(kind(0.d0)) :: grad_T - real(kind(0.d0)) :: tmp1, tmp2 real(kind(0.d0)) :: f_bpres_dot if (thermal == 3) then diff --git a/src/simulation/m_chemistry.fpp b/src/simulation/m_chemistry.fpp index 4465a8a4d..14b6c3b8e 100644 --- a/src/simulation/m_chemistry.fpp +++ b/src/simulation/m_chemistry.fpp @@ -55,8 +55,6 @@ contains type(vector_field), dimension(:), intent(IN) :: flux_n type(scalar_field), dimension(sys_size), intent(INOUT) :: rhs_vf - type(int_bounds_info) :: ix, iy, iz - integer :: x, y, z integer :: eqn @@ -103,21 +101,12 @@ contains subroutine s_compute_chemistry_reaction_flux(rhs_vf, q_cons_qp, q_prim_qp) type(scalar_field), dimension(sys_size), intent(INOUT) :: rhs_vf, q_cons_qp, q_prim_qp - - integer :: i - integer :: x, y, z integer :: eqn - real(kind(0d0)) :: T - integer :: o - real(kind(0d0)) :: dyn_pres - real(kind(0d0)) :: E - real(kind(0d0)) :: rho, omega_m real(kind(0d0)), dimension(num_species) :: Ys real(kind(0d0)), dimension(num_species) :: omega - real(kind(0d0)) :: cp_mix if (chemistry) then !$acc parallel loop collapse(3) gang vector default(present) & diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index dcba87d60..5c04fa28f 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -224,7 +224,6 @@ contains type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf integer, intent(IN) :: t_step - real(kind(0d0)), dimension(num_fluids) :: alpha_rho !< Cell-avg. partial density real(kind(0d0)) :: rho !< Cell-avg. density real(kind(0d0)), dimension(num_dims) :: vel !< Cell-avg. velocity real(kind(0d0)) :: vel_sum !< Cell-avg. velocity sum @@ -232,34 +231,17 @@ contains real(kind(0d0)), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction real(kind(0d0)) :: gamma !< Cell-avg. sp. heat ratio real(kind(0d0)) :: pi_inf !< Cell-avg. liquid stiffness function - real(kind(0d0)) :: qv !< Cell-avg. fluid reference energy real(kind(0d0)) :: c !< Cell-avg. sound speed - real(kind(0d0)) :: E !< Cell-avg. energy real(kind(0d0)) :: H !< Cell-avg. enthalpy real(kind(0d0)), dimension(2) :: Re !< Cell-avg. Reynolds numbers - - ! ICFL, VCFL, CCFL and Rc stability criteria extrema for the current - ! time-step and located on both the local (loc) and the global (glb) - ! computational domains - - real(kind(0d0)) :: blkmod1, blkmod2 !< - !! Fluid bulk modulus for Woods mixture sound speed - - integer :: i, j, k, l, q !< Generic loop iterators - - integer :: Nfq - real(kind(0d0)) :: fltr_dtheta !< - !! Modified dtheta accounting for Fourier filtering in azimuthal direction. + integer :: j, k, l ! Computing Stability Criteria at Current Time-step ================ - !$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho, vel, alpha, Re, fltr_dtheta, Nfq) + !$acc parallel loop collapse(3) gang vector default(present) private(vel, alpha, Re) do l = 0, p do k = 0, n do j = 0, m - call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l) - - ! Compute mixture sound speed call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0d0, c) if (viscous) then @@ -380,14 +362,9 @@ contains character(LEN=15) :: FMT - integer :: i, j, k, l, ii, r!< Generic loop iterators + integer :: i, j, k, l, r - real(kind(0d0)), dimension(nb) :: nRtmp !< Temporary bubble concentration - real(kind(0d0)) :: nbub, nR3, vftmp !< Temporary bubble number density real(kind(0d0)) :: gamma, lit_gamma, pi_inf, qv !< Temporary EOS params - real(kind(0d0)) :: rho !< Temporary density - real(kind(0d0)), dimension(2) :: Re !< Temporary Reynolds number - real(kind(0d0)) :: E_e !< Temp. elastic energy contribution ! Creating or overwriting the time-step root directory write (t_step_dir, '(A,I0,A,I0)') trim(case_dir)//'/p_all' @@ -950,21 +927,17 @@ contains real(kind(0d0)) :: int_pres real(kind(0d0)) :: max_pres real(kind(0d0)), dimension(2) :: Re - real(kind(0d0)) :: E_e real(kind(0d0)), dimension(6) :: tau_e real(kind(0d0)) :: G real(kind(0d0)) :: dyn_p, Temp - integer :: i, j, k, l, s, q, d !< Generic loop iterator + integer :: i, j, k, l, s, d !< Generic loop iterator real(kind(0d0)) :: nondim_time !< Non-dimensional time real(kind(0d0)) :: tmp !< !! Temporary variable to store quantity for mpi_allreduce - real(kind(0d0)) :: blkmod1, blkmod2 !< - !! Fluid bulk modulus for Woods mixture sound speed - integer :: npts !< Number of included integral points real(kind(0d0)) :: rad, thickness !< For integral quantities logical :: trigger !< For integral quantities @@ -978,7 +951,7 @@ contains if (t_step_old /= dflt_int) then nondim_time = real(t_step + t_step_old, kind(0d0))*dt else - nondim_time = real(t_step, kind(0d0))*dt !*1.d-5/10.0761131451d0 + nondim_time = real(t_step, kind(0d0))*dt end if end if @@ -1603,10 +1576,6 @@ contains !! other procedures that are necessary to setup the module. subroutine s_initialize_data_output_module - type(int_bounds_info) :: ix, iy, iz - - integer :: i !< Generic loop iterator - ! Allocating/initializing ICFL, VCFL, CCFL and Rc stability criteria @:ALLOCATE_GLOBAL(icfl_sf(0:m, 0:n, 0:p)) icfl_max = 0d0 @@ -1644,8 +1613,6 @@ contains !> Module deallocation and/or disassociation procedures subroutine s_finalize_data_output_module - integer :: i !< Generic loop iterator - ! Deallocating the ICFL, VCFL, CCFL, and Rc stability criteria @:DEALLOCATE_GLOBAL(icfl_sf) if (viscous) then diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index bc0552b1d..2552e37c2 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -871,10 +871,8 @@ contains logical :: beg_end_geq_0 - integer :: pack_offsets(1:3), unpack_offsets(1:3) integer :: pack_offset, unpack_offset real(kind(0d0)), pointer :: p_send, p_recv - integer, pointer, dimension(:) :: p_i_send, p_i_recv #ifdef MFC_MPI @@ -2136,7 +2134,6 @@ contains logical :: beg_end_geq_0 - integer :: pack_offsets(1:3), unpack_offsets(1:3) integer :: pack_offset, unpack_offset real(kind(0d0)), pointer :: p_send, p_recv diff --git a/src/simulation/m_qbmm.fpp b/src/simulation/m_qbmm.fpp index c478495a5..96e75a47d 100644 --- a/src/simulation/m_qbmm.fpp +++ b/src/simulation/m_qbmm.fpp @@ -702,7 +702,7 @@ contains real(kind(0.d0)), intent(in) :: pres, rho, c real(kind(0.d0)), dimension(nterms, 0:2, 0:2), intent(out) :: coeffs - integer :: i1, i2, q + integer :: i1, i2 coeffs = 0d0 @@ -776,7 +776,7 @@ contains real(kind(0.d0)), intent(inout) :: pres, rho, c real(kind(0.d0)), dimension(nterms, 0:2, 0:2), intent(out) :: coeffs - integer :: i1, i2, q + integer :: i1, i2 coeffs = 0d0 @@ -841,20 +841,19 @@ contains real(kind(0d0)), dimension(nmom) :: moms, msum real(kind(0d0)), dimension(nnode, nb) :: wght, abscX, abscY, wght_pb, wght_mv, wght_ht, ht - real(kind(0d0)), dimension(nterms, 0:2, 0:2) :: mom3d_terms, coeff - real(kind(0d0)) :: pres, rho, nbub, c, alf, R3, momsum, drdt, drdt2, chi_vw, x_vw, rho_mw, k_mw, T_bar, grad_T - real(kind(0d0)) :: start, finish + real(kind(0d0)), dimension(nterms, 0:2, 0:2) :: coeff + real(kind(0d0)) :: pres, rho, nbub, c, alf, momsum, drdt, drdt2, chi_vw, x_vw, rho_mw, k_mw, T_bar, grad_T real(kind(0d0)) :: n_tait, B_tait - integer :: j, k, l, q, r, s !< Loop variables integer :: id1, id2, id3 integer :: i1, i2 + integer :: j, q, r is1_qbmm = ix; is2_qbmm = iy; is3_qbmm = iz !$acc update device(is1_qbmm, is2_qbmm, is3_qbmm) - !$acc parallel loop collapse(3) gang vector default(present) private(moms, msum, wght, abscX, abscY, wght_pb, wght_mv, wght_ht, coeff, ht, r, q, n_tait, B_tait, pres, rho, nbub, c, alf, R3, momsum, drdt, drdt2, chi_vw, x_vw, rho_mw, k_mw, T_bar, grad_T) + !$acc parallel loop collapse(3) gang vector default(present) private(moms, msum, wght, abscX, abscY, wght_pb, wght_mv, wght_ht, coeff, ht, r, q, n_tait, B_tait, pres, rho, nbub, c, alf, momsum, drdt, drdt2, chi_vw, x_vw, rho_mw, k_mw, T_bar, grad_T) do id3 = is3_qbmm%beg, is3_qbmm%end do id2 = is2_qbmm%beg, is2_qbmm%end do id1 = is1_qbmm%beg, is1_qbmm%end diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index 084649837..37ab66476 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -699,30 +699,9 @@ contains integer, intent(in) :: t_step real(kind(0d0)), intent(inout) :: time_avg - real(kind(0d0)) :: t_start, t_finish - real(kind(0d0)) :: gp_sum - - real(kind(0d0)) :: top, bottom !< Numerator and denominator when evaluating flux limiter function - real(kind(0d0)), dimension(num_fluids) :: myalpha_rho, myalpha - - real(kind(0d0)) :: tmp1, tmp2, tmp3, tmp4, & - c_gas, c_liquid, & - Cpbw, Cpinf, Cpinf_dot, & - myH, myHdot, rddot, alf_gas - - real(kind(0d0)) :: n_tait, B_tait, angle, angle_z - - real(kind(0d0)), dimension(nb) :: Rtmp, Vtmp - real(kind(0d0)) :: myR, myV, alf, myP, myRho, R2Vav real(kind(0d0)), dimension(0:m, 0:n, 0:p) :: nbub - integer :: ndirs - - real(kind(0d0)) :: sound - real(kind(0d0)) :: start, finish - real(kind(0d0)) :: s2, const_sos, s1 - - integer :: i, c, j, k, l, q, ii, id !< Generic loop iterators - integer :: term_index + real(kind(0d0)) :: t_start, t_finish + integer :: i, j, k, l, id !< Generic loop iterators call nvtxStartRange("Compute_RHS") @@ -1593,7 +1572,7 @@ contains type(scalar_field), dimension(sys_size), intent(in) :: flux_src_n type(scalar_field), dimension(sys_size), intent(in) :: dq_prim_dx_vf, dq_prim_dy_vf, dq_prim_dz_vf - integer :: i, j, k, l, q + integer :: i, j, k, l if (idir == 1) then ! x-direction @@ -2089,7 +2068,6 @@ contains integer :: weno_dir !< Coordinate direction of the WENO reconstruction - integer :: i, j, k, l ! Reconstruction in s1-direction =================================== if (norm_dir == 1) then @@ -2213,7 +2191,7 @@ contains !> Module deallocation and/or disassociation procedures subroutine s_finalize_rhs_module - integer :: i, j, k, l !< Generic loop iterators + integer :: i, j, l do j = cont_idx%beg, cont_idx%end !$acc exit data detach(q_prim_qp%vf(j)%sf) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 99da073bb..f609c9593 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -314,7 +314,7 @@ contains real(kind(0d0)), dimension(num_species) :: Ys_L, Ys_R real(kind(0d0)), dimension(num_species) :: Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR real(kind(0d0)), dimension(num_species) :: Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2 - real(kind(0d0)) :: Cp_avg, Cv_avg, gamma_avggg, T_avg, eps, c_sum_Yi_Phi + real(kind(0d0)) :: Cp_avg, Cv_avg, T_avg, eps, c_sum_Yi_Phi real(kind(0d0)) :: T_L, T_R real(kind(0d0)) :: Y_L, Y_R real(kind(0d0)) :: MW_L, MW_R @@ -331,20 +331,15 @@ contains real(kind(0d0)), dimension(2) :: Re_L, Re_R real(kind(0d0)) :: rho_avg - real(kind(0d0)), dimension(num_dims) :: vel_avg real(kind(0d0)) :: H_avg real(kind(0d0)) :: gamma_avg real(kind(0d0)) :: c_avg real(kind(0d0)) :: s_L, s_R, s_M, s_P, s_S - real(kind(0d0)) :: xi_L, xi_R !< Left and right wave speeds functions real(kind(0d0)) :: xi_M, xi_P - real(kind(0d0)) :: nbub_L, nbub_R real(kind(0d0)) :: ptilde_L, ptilde_R real(kind(0d0)) :: vel_L_rms, vel_R_rms, vel_avg_rms - real(kind(0d0)) :: blkmod1, blkmod2 - real(kind(0d0)) :: rho_Star, E_Star, p_Star, p_K_Star real(kind(0d0)) :: Ms_L, Ms_R, pres_SL, pres_SR real(kind(0d0)) :: alpha_L_sum, alpha_R_sum @@ -374,7 +369,7 @@ contains if (norm_dir == ${NORM_DIR}$) then !$acc parallel loop collapse(3) gang vector default(present) & !$acc private(alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, & - !$acc alpha_R, vel_avg, tau_e_L, tau_e_R, G_L, G_R, Re_L, Re_R, & + !$acc alpha_R, tau_e_L, tau_e_R, G_L, G_R, Re_L, Re_R, & !$acc rho_avg, h_avg, gamma_avg, s_L, s_R, s_S, Ys_L, Ys_R, & !$acc Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, & !$acc Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2) @@ -947,7 +942,6 @@ contains real(kind(0d0)), dimension(2) :: Re_L, Re_R real(kind(0d0)) :: rho_avg - real(kind(0d0)), dimension(num_dims) :: vel_avg real(kind(0d0)) :: H_avg real(kind(0d0)) :: gamma_avg real(kind(0d0)) :: c_avg @@ -961,7 +955,6 @@ contains real(kind(0d0)), dimension(nb) :: V0_L, V0_R real(kind(0d0)), dimension(nb) :: P0_L, P0_R real(kind(0d0)), dimension(nb) :: pbw_L, pbw_R - real(kind(0d0)), dimension(nb, nmom) :: moms_L, moms_R real(kind(0d0)) :: ptilde_L, ptilde_R real(kind(0d0)) :: alpha_L_sum, alpha_R_sum, nbub_L_denom, nbub_R_denom @@ -972,10 +965,8 @@ contains real(kind(0d0)) :: vel_L_rms, vel_R_rms, vel_avg_rms real(kind(0d0)) :: vel_L_tmp, vel_R_tmp - real(kind(0d0)) :: blkmod1, blkmod2 real(kind(0d0)) :: rho_Star, E_Star, p_Star, p_K_Star real(kind(0d0)) :: pres_SL, pres_SR, Ms_L, Ms_R - real(kind(0d0)) :: start, finish real(kind(0d0)) :: zcoef, pcorr !< low Mach number correction integer :: i, j, k, l, q !< Generic loop iterators integer :: idx1, idxi @@ -1360,7 +1351,7 @@ contains end do elseif (model_eqns == 4) then !ME4 - !$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, vel_avg, & + !$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, & !$acc rho_avg, h_avg, gamma_avg, s_L, s_R, s_S, vel_avg_rms, nbub_L, nbub_R, ptilde_L, ptilde_R) do l = is3%beg, is3%end do k = is2%beg, is2%end diff --git a/src/simulation/m_sim_helpers.f90 b/src/simulation/m_sim_helpers.f90 index a4af91281..21dd369c3 100644 --- a/src/simulation/m_sim_helpers.f90 +++ b/src/simulation/m_sim_helpers.f90 @@ -86,7 +86,7 @@ end subroutine s_compute_enthalpy subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl_sf, vcfl_sf, Rc_sf) !$acc routine seq real(kind(0d0)), dimension(num_dims) :: vel - real(kind(0d0)) :: c, icfl_dt, vcfl_dt, rho + real(kind(0d0)) :: c, rho real(kind(0d0)), dimension(0:m, 0:n, 0:p) :: icfl_sf real(kind(0d0)), dimension(0:m, 0:n, 0:p), optional :: vcfl_sf, Rc_sf real(kind(0d0)) :: fltr_dtheta !< diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 633f0e8ce..0c27b1975 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -216,9 +216,6 @@ contains ! Logical used to check the existence of the current directory file logical :: file_exist - ! Generic loop iterators - integer :: i, j - ! Logistics ======================================================== file_path = trim(case_dir)//'/.' @@ -417,8 +414,6 @@ contains end if if (patch_ib(i)%c > 0) then - - print *, "HERE Np", Np allocate (airfoil_grid_u(1:Np)) allocate (airfoil_grid_l(1:Np)) @@ -1100,7 +1095,7 @@ contains real(kind(0d0)), intent(inout) :: start, finish integer, intent(inout) :: nt - integer :: i, j, k, l + integer :: i if (cfl_dt) then if (cfl_const_dt .and. t_step == 0) call s_compute_dt() diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index e9b260370..091bdbe2b 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -301,15 +301,9 @@ contains integer, intent(in) :: t_step real(kind(0d0)), intent(inout) :: time_avg - integer :: i, j, k, l, q!< Generic loop iterator - real(kind(0d0)) :: nR3bar - real(kind(0d0)) :: e_mix - - real(kind(0d0)) :: T - real(kind(0d0)), dimension(num_species) :: Ys + integer :: i, j, k, l, q !< Generic loop iterator ! Stage 1 of 1 ===================================================== - call nvtxStartRange("Time_Step") call s_compute_rhs(q_cons_ts(1)%vf, q_prim_vf, rhs_vf, pb_ts(1)%sf, rhs_pb, mv_ts(1)%sf, rhs_mv, t_step, time_avg) @@ -417,7 +411,6 @@ contains integer :: i, j, k, l, q!< Generic loop iterator real(kind(0d0)) :: start, finish - real(kind(0d0)) :: nR3bar ! Stage 1 of 2 ===================================================== @@ -599,9 +592,7 @@ contains real(kind(0d0)), intent(INOUT) :: time_avg integer :: i, j, k, l, q !< Generic loop iterator - real(kind(0d0)) :: ts_error, denom, error_fraction, time_step_factor !< Generic loop iterator real(kind(0d0)) :: start, finish - real(kind(0d0)) :: nR3bar ! Stage 1 of 3 ===================================================== @@ -862,7 +853,6 @@ contains integer, intent(in) :: t_step real(kind(0d0)), intent(inout) :: time_avg - integer :: i, j, k, l !< Generic loop iterator real(kind(0d0)) :: start, finish call cpu_time(start) @@ -896,8 +886,6 @@ contains type(vector_field) :: gm_alpha_qp - integer :: i, j, k, l, q !< Generic loop iterator - call s_convert_conservative_to_primitive_variables( & q_cons_ts(1)%vf, & q_prim_vf, & @@ -924,7 +912,7 @@ contains real(kind(0d0)), dimension(2) :: Re !< Cell-avg. Reynolds numbers type(vector_field) :: gm_alpha_qp real(kind(0d0)) :: dt_local - integer :: i, j, k, l, q !< Generic loop iterators + integer :: j, k, l !< Generic loop iterators call s_convert_conservative_to_primitive_variables( & q_cons_ts(1)%vf, & diff --git a/src/simulation/m_viscous.fpp b/src/simulation/m_viscous.fpp index 3e2c3bf70..32942ed1c 100644 --- a/src/simulation/m_viscous.fpp +++ b/src/simulation/m_viscous.fpp @@ -1467,8 +1467,6 @@ contains subroutine s_finalize_viscous_module() - integer :: i - @:DEALLOCATE_GLOBAL(Res_viscous) end subroutine s_finalize_viscous_module diff --git a/src/simulation/m_weno.fpp b/src/simulation/m_weno.fpp index f0b5490b4..655e0c6a5 100644 --- a/src/simulation/m_weno.fpp +++ b/src/simulation/m_weno.fpp @@ -160,7 +160,6 @@ contains !! other procedures that are necessary to setup the module. subroutine s_initialize_weno_module - integer :: i, j if (weno_order == 1) return ! Allocating/Computing WENO Coefficients in x-direction ============ @@ -552,7 +551,6 @@ contains real(kind(0d0)), dimension(0:weno_num_stencils) :: delta real(kind(0d0)), dimension(-3:3) :: v ! temporary field value array for clarity (WENO7 only) real(kind(0d0)) :: tau - real(kind(0d0)), pointer :: beta_p(:) integer :: i, j, k, l @@ -974,7 +972,7 @@ contains integer, intent(IN) :: norm_dir integer, intent(IN) :: weno_dir - integer :: i, j, k, l, q !< Generic loop iterators + integer :: j, k, l, q ! Determining the number of cell-average variables which will be ! WENO-reconstructed and mapping their indical bounds in the x-, @@ -1253,8 +1251,6 @@ contains !> Module deallocation and/or disassociation procedures subroutine s_finalize_weno_module() - integer :: i, j - if (weno_order == 1) return ! Deallocating the WENO-stencil of the WENO-reconstructed variables From 6f09c065f000d95111760dada71322ecd287d5cb Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 8 Nov 2024 09:46:43 -0500 Subject: [PATCH 8/9] Move codecov setting (#710) --- .github/{workflows => }/codecov.yml | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename .github/{workflows => }/codecov.yml (100%) diff --git a/.github/workflows/codecov.yml b/.github/codecov.yml similarity index 100% rename from .github/workflows/codecov.yml rename to .github/codecov.yml From 9700eb5ab2a5b7906db8bcfc0f08380726e9485d Mon Sep 17 00:00:00 2001 From: Max Hawkins <37495064+max-Hawkins@users.noreply.github.com> Date: Fri, 8 Nov 2024 11:00:31 -0500 Subject: [PATCH 9/9] Add Phoenix Quadro RTX6000 CUDA CC 75 (#711) --- toolchain/modules | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/toolchain/modules b/toolchain/modules index 2d0ebe925..ca061b025 100644 --- a/toolchain/modules +++ b/toolchain/modules @@ -45,7 +45,7 @@ p GT Phoenix p-all python/3.10.10 p-cpu gcc/12.3.0 openmpi/4.1.5 p-gpu nvhpc/24.5 hpcx/2.19-cuda cuda/12.1.1 -p-gpu MFC_CUDA_CC=70,80,89,90 NVHPC_CUDA_HOME=$CUDA_HOME CC=nvc CXX=nvc++ FC=nvfortran +p-gpu MFC_CUDA_CC=70,75,80,89,90 NVHPC_CUDA_HOME=$CUDA_HOME CC=nvc CXX=nvc++ FC=nvfortran f OLCF Frontier f-all cce/18.0.0 cpe/24.07 rocm/6.1.3 cray-mpich/8.1.28