diff --git a/CHANGELOG.md b/CHANGELOG.md index 033e17f1807d..d416562464bf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,11 @@ # Changelog ## Develop +- *BREAKING* Modularized entropy viscosity in the compressible Navier-Stokes + solver. The obsolete `case.numerics.c_avisc_low` and + `case.numerics.c_avisc_entropy` options are now rejected. Configure an + `artificial_viscosity_model` simulation component with those coefficients + and enable its field through the fluid `viscous_regularization` object. - *BREAKING*, normal_outflow conditions now require specifying `value`, which is used to set the value of the tangential components of velocity. @@ -153,6 +158,9 @@ - Added `host_array_t` and `device_array_t` temporary array types and support for requesting these through `scratch_registry_t`. - Added the `cai_sagaut_model_ii` wall model with CPU, CUDA, HIP, and OpenCL. +- Modularized `entropy_viscosity` into two steps: 1. Calculate the + artificial viscosity field in a simcomp; 2. Apply the artificial + viscosity though a viscous regularization object. - Added the `create_periodic_zones` contrib utility for converting pairs of labeled zones in an existing `.nmsh` mesh into periodic zones. backends. This model is based on the work of [Cai and Sagaut (PoF, diff --git a/doc/pages/user-guide/case-file.md b/doc/pages/user-guide/case-file.md index e2c696377d44..57d11de2fde2 100644 --- a/doc/pages/user-guide/case-file.md +++ b/doc/pages/user-guide/case-file.md @@ -373,13 +373,6 @@ The compressible solver accepts the following parameters: | `mu` | Constant physical dynamic viscosity | Non-negative real | `0.0` | | `kappa` | Constant physical thermal conductivity | Non-negative real | `0.0` | -Additional numerics parameters specific to compressible flows: - -| Name | Description | Admissible values | Default value | -| ----------------- | -------------------------------------------------- | ----------------- | ------------- | -| `c_avisc_low` | Coefficient for low-order artificial viscosity | Positive reals | `0.5` | -| `c_avisc_entropy` | Coefficient for entropy-based artificial viscosity | Positive reals | `1.0` | - The compressible solver uses variable time-stepping controlled by the CFL number. Set `variable_timestep` to `true` and specify `target_cfl` in the time control object. @@ -403,6 +396,10 @@ Example configuration: "fluid": { "scheme": "compressible", "gamma": 1.4, + "viscous_regularization":{ + "type": "artificial_viscosity", + "reg_coeff_name": "entropy_viscosity" + }, "initial_condition": { "type": "user" }, @@ -426,12 +423,13 @@ Example configuration: "output_control": "nsamples", "output_value": 20 }, - "numerics": { - "time_order": 3, - "polynomial_order": 5, - "c_avisc_low": 0.5, - "c_avisc_entropy": 0.5 - } + "simulation_components": [ + { + "type": "artificial_viscosity_model", + "model": "entropy_viscosity", + "field_name": "entropy_viscosity" + } + ] } ~~~~~~~~~~~~~~~ @@ -2020,7 +2018,6 @@ concisely directly in the table. | `shear_stress.value` | The shear stress vector value for `sh` boundaries | Vector of 3 reals | `[0, 0, 0]` | | `wall_modelling.type` | The wall model type for `wm` boundaries. See documentation for additional config parameters. | `rough_log_law`, `spalding` | - | | `source_terms` | Array of JSON objects, defining additional source terms. | See list of source terms above | - | -| `gradient_jump_penalty` | Array of JSON objects, defining additional gradient jump penalty. | See list of gradient jump penalty above | - | | `boundary_types` | Boundary types/conditions labels. | Array of strings | - | | `velocity_solver.type` | Linear solver for the momentum equation. | `cg`, `pipecg`, `bicgstab`, `cacg`, `gmres` | - | | `velocity_solver.preconditioner.type` | Linear solver preconditioner for the momentum equation. | `ident`, `hsmg`, `jacobi` | - | @@ -2237,3 +2234,43 @@ currently supports 50 regions, with id 1..25 being reserved for internal use. | ---------------- | ----------------------------------------------------------- | ----------------- | ------------- | | `enabled` | Whether to enable gathering of runtime statistics | `true` or `false` | `false` | | `output_profile` | Whether to output all gathered profiling data as a CSV file | `true` or `false` | `false` | + +## Viscous regularization {#case-file_viscous-regularization} + +Users can use the viscous regularization object to enhance the smoothness or +the numerical stability of the solution via a diffusion term. For example, +artificial viscosity can be set up by this object to perform shock capturing. +Note that some regularization techniques do not have a diffusive mathematical +form, for example gradient jump penalty and high-pass filter relaxation terms, +and they are included in the source terms instead. The viscous regularization +can be set up with the following options: + +* `type`, the viscous regularization type. + - `artificial_viscosity`, the standard second-order diffusion term, + $\frac{\partial}{\partial x}\left(\mu_\mathrm{artificial}\frac{\partial u}{\partial x}\right)$. +* `reg_coeff_name`, name of the $\mu_\mathrm{artificial}$ field, usually +computed by a simulation component. + +Viscous regularization is currently supported by the compressible fluid +solver. Artificial viscosity requires both the regularization object that +consumes the coefficient and a simulation component that computes it. The +`reg_coeff_name` and `field_name` values must match: + +~~~~~~~~~~~~~~~{.json} +{ + "fluid": { + "scheme": "compressible", + "viscous_regularization": { + "type": "artificial_viscosity", + "reg_coeff_name": "entropy_viscosity" + } + }, + "simulation_components": [ + { + "type": "artificial_viscosity_model", + "model": "entropy_viscosity", + "field_name": "entropy_viscosity" + } + ] +} +~~~~~~~~~~~~~~~ diff --git a/doc/pages/user-guide/simcomps.md b/doc/pages/user-guide/simcomps.md index cd7dadd0ddbe..ad1ab2ec1f7c 100644 --- a/doc/pages/user-guide/simcomps.md +++ b/doc/pages/user-guide/simcomps.md @@ -65,6 +65,8 @@ in Neko. The list will be updated as new simcomps are added. - Streaming of data for in-situ field manipulation \ref simcomp_data_streamer - Sub-sampling of fields by changing polynomial order and masking by point zones \ref simcomp_field_subsampler +- Computation of the viscous coefficient in artificial viscosity method (AVM) + via a particular model \ref simcomp_avm ## Controlling execution and file output Each simulation component is, by default, executed once per time step to perform @@ -764,15 +766,15 @@ keywords: the maximum spacing of the GLL nodes within the element. The `les_delta` field is added to the registry and written to the .fld files. - `nut_field`: The name of the SGS eddy viscosity field added to the registry. - Defaults to `nut`. This allows to have two different SGS models active, saved - to different fields. For example, one for the scalar and one to the fluid. -- `extrapolation`: Whether or not extrapolate the velocity to + Defaults to `nut`. This allows two different SGS models to be active and saved + to different fields, for example one for the scalar and one for the fluid. +- `extrapolation`: Whether or not to extrapolate the velocity to compute the eddy viscosity. - - `true`: extrapolate the velocity as the same order as + - `true`: extrapolate the velocity to the same order as the time scheme. - - `false`: the default option, disable the extrapolation. + - `false`: the default option, disables extrapolation. In this case, the estimation of the eddy viscosity is of first order, while - circumvent the risk of unstable extrapolation. + circumventing the risk of unstable extrapolation. ~~~~~~~~~~~~~~~{.json} { @@ -932,4 +934,34 @@ keywords used by the latter can also be specified, with the exception of "compute_control": "tsteps", "compute_value": 10 } - ~~~~~~~~~~~~~~~ \ No newline at end of file + ~~~~~~~~~~~~~~~ + +### artificial_viscosity_model {#simcomp_avm} +Computes an artificial viscosity field using a particular model. **Note:** The simcomp +*only* computes the artificial viscosity field. You have to select the +corresponding `reg_coeff_name` in the compressible fluid solver's +`viscous_regularization` JSON object to actually enable artificial viscosity; +see @ref case-file_viscous-regularization. The simcomp is controlled by the +following keywords: + +- `model`: Selects the artificial viscosity model. Currently available models are: + - `entropy_viscosity`: An entropy-based viscosity aimed at shock capturing; + see Guermond et al. (https://doi.org/10.1016/j.jcp.2010.11.043). + - `c_avisc_low`: Coefficient for low-order artificial viscosity, defaults + to 0.5 + - `c_avisc_entropy`: Coefficient for entropy-based artificial viscosity, + defaults to 1.0 + +- `field_name`: The name of the artificial viscosity field added to the registry. + Defaults to the model name. + It must match the `reg_coeff_name` selected by the consuming regularization + object. + + ~~~~~~~~~~~~~~~{.json} + { + "type": "artificial_viscosity_model", + "name": "artificial_viscosity_model", + "model": "entropy_viscosity", + "field_name": "entropy_viscosity" + } + ~~~~~~~~~~~~~~~ diff --git a/doc/schemas/numerics.schema.json b/doc/schemas/numerics.schema.json index 5d7689d268f5..aaff5f506355 100644 --- a/doc/schemas/numerics.schema.json +++ b/doc/schemas/numerics.schema.json @@ -43,12 +43,6 @@ }, "oifs_target_cfl": { "$ref": "urn:neko:schema:common#/$defs/positiveNumberOrRef" - }, - "c_avisc_low": { - "$ref": "urn:neko:schema:common#/$defs/positiveNumberOrRef" - }, - "c_avisc_entropy": { - "$ref": "urn:neko:schema:common#/$defs/positiveNumberOrRef" } }, "additionalProperties": false diff --git a/doc/schemas/simcomps/artificial_viscosity_model.schema.json b/doc/schemas/simcomps/artificial_viscosity_model.schema.json new file mode 100644 index 000000000000..c51992823e62 --- /dev/null +++ b/doc/schemas/simcomps/artificial_viscosity_model.schema.json @@ -0,0 +1,40 @@ +{ + "$schema": "https://json-schema.org/draft/2020-12/schema", + "$id": "urn:neko:schema:simcomps:artificial_viscosity_model", + "$defs": { + "simulationComponent": { + "allOf": [ + { + "$ref": "urn:neko:schema:simulation-components#/$defs/simcompBase" + }, + { + "$ref": "urn:neko:schema:simulation-components#/$defs/simcompFieldWriterOptions" + }, + { + "type": "object", + "required": [ + "model" + ], + "properties": { + "type": { + "const": "artificial_viscosity_model" + }, + "model": { + "const": "entropy_viscosity" + }, + "field_name": { + "$ref": "urn:neko:schema:simulation-components#/$defs/nonEmptyString" + }, + "c_avisc_low": { + "$ref": "urn:neko:schema:simulation-components#/$defs/numberOrRef" + }, + "c_avisc_entropy": { + "$ref": "urn:neko:schema:simulation-components#/$defs/numberOrRef" + } + }, + "additionalProperties": true + } + ] + } + } +} diff --git a/doc/schemas/simulation-components.schema.json b/doc/schemas/simulation-components.schema.json index 8dab25d70454..0fef4e2f4196 100644 --- a/doc/schemas/simulation-components.schema.json +++ b/doc/schemas/simulation-components.schema.json @@ -607,6 +607,9 @@ { "$ref": "urn:neko:schema:simcomps:field_subsampler#/$defs/simulationComponent" }, + { + "$ref": "urn:neko:schema:simcomps:artificial_viscosity_model#/$defs/simulationComponent" + }, { "$ref": "urn:neko:schema:simcomps:boundary_data_writer#/$defs/simulationComponent" }, @@ -616,4 +619,4 @@ ] } } -} \ No newline at end of file +} diff --git a/examples/cns_2d_shock_tube/2d_shock_tube.case b/examples/cns_2d_shock_tube/2d_shock_tube.case index 3bbdc6e0406a..238c9a9308ea 100644 --- a/examples/cns_2d_shock_tube/2d_shock_tube.case +++ b/examples/cns_2d_shock_tube/2d_shock_tube.case @@ -13,13 +13,15 @@ }, "numerics": { "time_order": 3, - "polynomial_order": 3, - "c_avisc_low": 0.5, - "c_avisc_entropy": 0.5 + "polynomial_order": 3 }, "fluid": { "scheme": "compressible", "gamma": 1.4, + "viscous_regularization": { + "type": "artificial_viscosity", + "reg_coeff_name": "entropy_viscosity" + }, "Re": 1000.0, "Pr": 0.73, "initial_condition": { @@ -39,6 +41,12 @@ "output_value": 100 }, "simulation_components": [ + { + "type": "artificial_viscosity_model", + "model": "entropy_viscosity", + "c_avisc_low": 0.5, + "c_avisc_entropy": 0.5 + }, { "type": "field_writer", "fields": ["fluid_rho", "S", "artificial_visc"], diff --git a/examples/euler_1d_sod/sod.case b/examples/euler_1d_sod/sod.case index 49e7c993fa12..338ae6ef7605 100755 --- a/examples/euler_1d_sod/sod.case +++ b/examples/euler_1d_sod/sod.case @@ -13,13 +13,15 @@ }, "numerics": { "time_order": 4, - "polynomial_order": 3, - "c_avisc_low": 0.5, - "c_avisc_entropy": 1.0 + "polynomial_order": 3 }, "fluid": { "scheme": "compressible", "gamma": 1.4, + "viscous_regularization":{ + "type": "artificial_viscosity", + "reg_coeff_name": "entropy_viscosity" + }, "initial_condition": { "type": "user" }, @@ -33,6 +35,10 @@ "output_value": 100 }, "simulation_components": [ + { + "type": "artificial_viscosity_model", + "model": "entropy_viscosity" + }, { "type": "field_writer", "fields": ["fluid_rho", "S", "artificial_visc"], // s2, s3, s4 in Paraview diff --git a/examples/euler_2d_cylinder/euler_2d_cylinder.case b/examples/euler_2d_cylinder/euler_2d_cylinder.case index 5eb52be08f68..ddd1cd2649c1 100644 --- a/examples/euler_2d_cylinder/euler_2d_cylinder.case +++ b/examples/euler_2d_cylinder/euler_2d_cylinder.case @@ -14,14 +14,16 @@ "numerics": { "time_order": 3, - "polynomial_order": 5, - "c_avisc_low": 0.5, - "c_avisc_entropy": 1.0 + "polynomial_order": 5 }, "fluid": { "scheme": "compressible", "gamma": 1.4, + "viscous_regularization":{ + "type": "artificial_viscosity", + "reg_coeff_name": "entropy_viscosity" + }, "initial_condition": { "type": "user" @@ -58,6 +60,10 @@ }, "simulation_components": [ + { + "type": "artificial_viscosity_model", + "model": "entropy_viscosity" + }, { "type": "field_writer", "fields": ["fluid_rho", "S", "artificial_visc"], // s2, s3, s4 in Paraview diff --git a/examples/euler_2d_forward_facing_step/step.case b/examples/euler_2d_forward_facing_step/step.case index e07f15bade65..e7989f4cd317 100644 --- a/examples/euler_2d_forward_facing_step/step.case +++ b/examples/euler_2d_forward_facing_step/step.case @@ -14,14 +14,16 @@ "numerics": { "time_order": 3, - "polynomial_order": 5, - "c_avisc_low": 0.5, - "c_avisc_entropy": 0.5 + "polynomial_order": 5 }, "fluid": { "scheme": "compressible", "gamma": 1.4, + "viscous_regularization":{ + "type": "artificial_viscosity", + "reg_coeff_name": "entropy_viscosity" + }, "initial_condition": { "type": "user" @@ -54,6 +56,10 @@ }, "simulation_components": [ + { + "type": "artificial_viscosity_model", + "model": "entropy_viscosity" + }, { "type": "field_writer", "fields": ["fluid_rho", "S", "artificial_visc", "max_wave_speed"], // s2, s3, s4, s5 in Paraview diff --git a/examples/euler_2d_smooth/euler_2d_smooth.case b/examples/euler_2d_smooth/euler_2d_smooth.case index 311611add25d..e1d053bc480f 100755 --- a/examples/euler_2d_smooth/euler_2d_smooth.case +++ b/examples/euler_2d_smooth/euler_2d_smooth.case @@ -13,13 +13,15 @@ }, "numerics": { "time_order": 3, - "polynomial_order": 3, - "c_avisc_low": 0.5, - "c_avisc_entropy": 1.0 + "polynomial_order": 3 }, "fluid": { "scheme": "compressible", "gamma": 1.4, + "viscous_regularization":{ + "type": "artificial_viscosity", + "reg_coeff_name": "entropy_viscosity" + }, "initial_condition": { "type": "user" }, @@ -27,6 +29,10 @@ "output_value": 100 }, "simulation_components": [ + { + "type": "artificial_viscosity_model", + "model": "entropy_viscosity" + }, { "type": "field_writer", "fields": ["fluid_rho", "S", "artificial_visc"], // s2, s3, s4 in Paraview diff --git a/examples/euler_tgv/euler_tgv.case b/examples/euler_tgv/euler_tgv.case index 92ba34c6eb16..681d3055e58a 100755 --- a/examples/euler_tgv/euler_tgv.case +++ b/examples/euler_tgv/euler_tgv.case @@ -13,13 +13,15 @@ }, "numerics": { "time_order": 3, - "polynomial_order": 5, - "c_avisc_low": 0.5, - "c_avisc_entropy": 1.0 + "polynomial_order": 5 }, "fluid": { "scheme": "compressible", "gamma": 1.4, + "viscous_regularization":{ + "type": "artificial_viscosity", + "reg_coeff_name": "entropy_viscosity" + }, "initial_condition": { "type": "user" }, @@ -27,6 +29,10 @@ "output_value": 100 }, "simulation_components": [ + { + "type": "artificial_viscosity_model", + "model": "entropy_viscosity" + }, { "type": "field_writer", "fields": ["fluid_rho", "S", "artificial_visc"], // s2, s3, s4 in Paraview diff --git a/src/.depends b/src/.depends index 2de49ba3b839..b22ffc1f0e6d 100644 --- a/src/.depends +++ b/src/.depends @@ -285,12 +285,7 @@ fluid/stress_formulation/bcknd/device/pnpn_res_stress_device.lo : fluid/stress_f fluid/fluid_scheme_base.lo : fluid/fluid_scheme_base.f90 common/time_state.lo field/field_list.lo bc/bc_list.lo common/utils.lo common/user_intf.lo common/time_step_controller.lo time_schemes/time_scheme_controller.lo sem/space.lo mesh/mesh.lo common/checkpoint.lo config/num_types.lo common/log.lo gs/gather_scatter.lo field/field_series.lo field/field.lo sem/dofmap.lo bc/dirichlet.lo sem/coef.lo fluid/fluid_scheme_incompressible.lo : fluid/fluid_scheme_incompressible.f90 common/time_state.lo device/device.lo math/field_math.lo common/utils.lo common/user_intf.lo registries/scratch_registry.lo common/json_utils.lo registries/registry.lo common/log.lo math/operators.lo math/math.lo mesh/mesh.lo bc/bc_list.lo bc/bc.lo fluid/fluid_stats.lo krylov/precon.lo multigrid/phmg.lo krylov/pc_hsmg.lo krylov/bcknd/device/pc_jacobi_device.lo krylov/bcknd/sx/pc_jacobi_sx.lo krylov/bcknd/cpu/pc_jacobi.lo bc/dirichlet.lo sem/coef.lo krylov/krylov.lo sem/dofmap.lo sem/space.lo field/field.lo fluid/fluid_source_term.lo config/num_types.lo common/checkpoint.lo config/neko_config.lo gs/gather_scatter.lo fluid/fluid_scheme_base.lo fluid/fluid_scheme_compressible.lo : fluid/fluid_scheme_compressible.f90 math/math.lo common/log.lo common/time_state.lo config/neko_config.lo fluid/bcknd/device/compressible_ops_device.lo fluid/bcknd/cpu/compressible_ops_cpu.lo device/device.lo math/operators.lo common/json_utils.lo common/user_intf.lo sem/space.lo registries/scratch_registry.lo mesh/mesh.lo config/num_types.lo fluid/fluid_scheme_base.lo bc/bc.lo registries/registry.lo math/field_math.lo field/field.lo -fluid/regularization.lo : fluid/regularization.f90 common/time_state.lo sem/dofmap.lo sem/coef.lo field/field.lo config/num_types.lo -fluid/regularization_fctry.lo : fluid/regularization_fctry.f90 common/utils.lo fluid/regularization/entropy_viscosity.lo fluid/regularization.lo -fluid/regularization/entropy_viscosity.lo : fluid/regularization/entropy_viscosity.f90 fluid/regularization/bcknd/device/entropy_viscosity_device.lo fluid/regularization/bcknd/cpu/entropy_viscosity_cpu.lo math/bcknd/device/device_math.lo device/device.lo config/neko_config.lo gs/gs_ops.lo gs/gather_scatter.lo sem/space.lo mesh/mesh.lo registries/scratch_registry.lo math/math.lo math/operators.lo time_schemes/bdf_time_scheme.lo common/time_state.lo sem/dofmap.lo sem/coef.lo field/field_series.lo math/field_math.lo field/field.lo common/json_utils.lo fluid/regularization.lo config/num_types.lo -fluid/regularization/bcknd/cpu/entropy_viscosity_cpu.lo : fluid/regularization/bcknd/cpu/entropy_viscosity_cpu.f90 config/num_types.lo -fluid/regularization/bcknd/device/entropy_viscosity_device.lo : fluid/regularization/bcknd/device/entropy_viscosity_device.F90 common/utils.lo config/num_types.lo -fluid/fluid_scheme_compressible_ns.lo : fluid/fluid_scheme_compressible_ns.f90 fluid/regularization/entropy_viscosity.lo bc/vector_bc_projector.lo fluid/regularization.lo fluid/bcknd/device/compressible_ops_device.lo fluid/bcknd/cpu/compressible_ops_cpu.lo common/time_state.lo common/log.lo common/utils.lo bc/bc.lo bc/bc_list.lo time_schemes/runge_kutta_scheme.lo config/neko_config.lo fluid/compressible_res.lo sem/coef.lo math/ax.lo common/time_step_controller.lo common/user_intf.lo common/profiler.lo common/json_utils.lo common/checkpoint.lo mesh/mesh.lo config/num_types.lo gs/gather_scatter.lo gs/gs_ops.lo registries/scratch_registry.lo fluid/fluid_scheme_compressible.lo field/field.lo math/bcknd/device/device_math.lo math/math.lo time_schemes/time_scheme_controller.lo time_schemes/bdf_time_scheme.lo field/field_series.lo math/operators.lo device/device.lo fluid/advection.lo comm/comm.lo +fluid/fluid_scheme_compressible_ns.lo : fluid/fluid_scheme_compressible_ns.f90 bc/vector_bc_projector.lo viscous_regularization/viscous_regularization.lo fluid/bcknd/device/compressible_ops_device.lo fluid/bcknd/cpu/compressible_ops_cpu.lo common/time_state.lo common/log.lo common/utils.lo bc/bc.lo bc/bc_list.lo time_schemes/runge_kutta_scheme.lo config/neko_config.lo fluid/compressible_res.lo sem/coef.lo math/ax.lo common/time_step_controller.lo common/user_intf.lo common/profiler.lo common/json_utils.lo common/checkpoint.lo mesh/mesh.lo config/num_types.lo gs/gather_scatter.lo gs/gs_ops.lo registries/scratch_registry.lo fluid/fluid_scheme_compressible.lo field/field.lo time_schemes/time_scheme_controller.lo time_schemes/bdf_time_scheme.lo field/field_series.lo math/operators.lo device/device.lo fluid/advection.lo comm/comm.lo fluid/fluid_aux.lo : fluid/fluid_aux.f90 common/time_state.lo common/utils.lo common/log.lo krylov/krylov.lo fluid/fluid_pnpn.lo : fluid/fluid_pnpn.f90 bc/no_slip.lo bc/dong_outflow.lo bc/field_dirichlet_vector.lo bc/blasius.lo bc/field_dirichlet.lo bc/inflow.lo field/field_series.lo ale/ale_manager.lo comm/comm.lo common/time_state.lo math/bcknd/device/opr_device.lo math/operators.lo io/file.lo bc/vector_bc_projector.lo bc/scalar_bc_projector.lo bc/mixed_bc.lo bc/bc.lo math/field_math.lo common/utils.lo bc/zero_dirichlet.lo math/mathops.lo config/neko_config.lo gs/gs_ops.lo common/time_step_controller.lo common/user_intf.lo mesh/mesh.lo common/checkpoint.lo bc/symmetry.lo bc/symmetry_aligned.lo bc/non_normal.lo bc/non_normal_aligned.lo bc/facet_normal.lo bc/wall_model_bc.lo bc/shear_stress.lo bc/dirichlet.lo field/field.lo math/ax.lo common/json_utils.lo common/profiler.lo fluid/advection.lo device/device.lo common/projection_vel.lo common/projection.lo fluid/fluid_aux.lo math/bcknd/device/device_mathops.lo registries/scratch_registry.lo fluid/fluid_scheme_incompressible.lo fluid/fluid_volflow.lo common/rhs_maker.lo fluid/pnpn_res.lo krylov/krylov.lo config/num_types.lo common/log.lo registries/registry.lo sem/coef.lo fluid/fluid_base_fctry.lo : fluid/fluid_base_fctry.f90 common/utils.lo fluid/fluid_scheme_compressible_ns.lo fluid/fluid_pnpn.lo fluid/fluid_scheme_base.lo @@ -400,7 +395,8 @@ simulation_components/weak_gradient_simcomp.lo : simulation_components/weak_grad simulation_components/gradient_simcomp.lo : simulation_components/gradient_simcomp.f90 common/utils.lo common/time_based_controller.lo simulation_components/field_writer.lo common/json_utils.lo case.lo common/time_state.lo math/operators.lo field/field.lo registries/registry.lo simulation_components/simulation_component.lo config/num_types.lo simulation_components/derivative_simcomp.lo : simulation_components/derivative_simcomp.f90 common/time_based_controller.lo common/utils.lo simulation_components/field_writer.lo common/json_utils.lo case.lo math/operators.lo field/field.lo registries/registry.lo common/time_state.lo simulation_components/simulation_component.lo config/num_types.lo simulation_components/les_simcomp.lo : simulation_components/les_simcomp.f90 common/utils.lo simulation_components/field_writer.lo common/json_utils.lo les/les_model.lo common/time_state.lo case.lo simulation_components/simulation_component.lo config/num_types.lo -simulation_components/simulation_component_fctry.lo : simulation_components/simulation_component_fctry.f90 simulation_components/boundary_data_writer_simcomp.lo simulation_components/wall_shear_stress_simcomp.lo simulation_components/field_subsampler.lo simulation_components/data_streamer_simcomp.lo simulation_components/spectral_error.lo simulation_components/derivative_simcomp.lo simulation_components/divergence_simcomp.lo simulation_components/gradient_simcomp.lo simulation_components/weak_gradient_simcomp.lo simulation_components/curl_simcomp.lo simulation_components/field_writer.lo common/utils.lo simulation_components/les_simcomp.lo simulation_components/probes.lo simulation_components/lambda2.lo simulation_components/user_stats.lo simulation_components/spatial_average.lo simulation_components/scalar_sgs_stats_simcomp.lo simulation_components/scalar_stats_simcomp.lo simulation_components/lpt_simcomp.lo simulation_components/fluid_sgs_stats_simcomp.lo simulation_components/fluid_stats_simcomp.lo simulation_components/force_torque.lo simulation_components/boundary_flux.lo simulation_components/boundary_operation.lo simulation_components/simulation_component.lo +simulation_components/simulation_component_fctry.lo : simulation_components/simulation_component_fctry.f90 simulation_components/boundary_data_writer_simcomp.lo simulation_components/wall_shear_stress_simcomp.lo simulation_components/avm_simcomp.lo simulation_components/field_subsampler.lo simulation_components/data_streamer_simcomp.lo simulation_components/spectral_error.lo simulation_components/derivative_simcomp.lo simulation_components/divergence_simcomp.lo simulation_components/gradient_simcomp.lo simulation_components/weak_gradient_simcomp.lo simulation_components/curl_simcomp.lo simulation_components/field_writer.lo common/utils.lo simulation_components/les_simcomp.lo simulation_components/probes.lo simulation_components/lambda2.lo simulation_components/user_stats.lo simulation_components/spatial_average.lo simulation_components/scalar_sgs_stats_simcomp.lo simulation_components/scalar_stats_simcomp.lo simulation_components/lpt_simcomp.lo simulation_components/fluid_sgs_stats_simcomp.lo simulation_components/fluid_stats_simcomp.lo simulation_components/force_torque.lo simulation_components/boundary_flux.lo simulation_components/boundary_operation.lo simulation_components/simulation_component.lo +simulation_components/avm_simcomp.lo : simulation_components/avm_simcomp.f90 common/utils.lo simulation_components/field_writer.lo common/json_utils.lo viscous_regularization/avm/avm.lo common/time_state.lo case.lo simulation_components/simulation_component.lo config/num_types.lo source_terms/source_term.lo : source_terms/source_term.f90 common/time_state.lo field/field_list.lo sem/coef.lo config/num_types.lo source_terms/user_source_term.lo : source_terms/user_source_term.f90 common/time_state.lo common/user_intf.lo sem/dofmap.lo math/field_math.lo sem/coef.lo field/field_list.lo source_terms/source_term.lo common/utils.lo config/num_types.lo config/neko_config.lo source_terms/translation_source_term.lo : source_terms/translation_source_term.f90 registries/scratch_registry.lo common/time_state.lo registries/registry.lo device/device.lo math/operators.lo math/field_math.lo field/field.lo common/utils.lo config/neko_config.lo sem/coef.lo source_terms/source_term.lo common/json_utils.lo field/field_list.lo config/num_types.lo @@ -485,4 +481,12 @@ wall_models/richardson.lo : wall_models/richardson.f90 math/bcknd/device/device_ wall_models/bcknd/cpu/richardson_cpu.lo : wall_models/bcknd/cpu/richardson_cpu.f90 math/math.lo common/log.lo common/utils.lo config/num_types.lo wall_models/bcknd/device/richardson_device.lo : wall_models/bcknd/device/richardson_device.F90 common/utils.lo config/num_types.lo wall_models/wall_model_fctry.lo : wall_models/wall_model_fctry.f90 common/utils.lo wall_models/richardson.lo wall_models/most.lo wall_models/rough_log_law.lo wall_models/cai_sagaut_model_ii.lo wall_models/spalding.lo wall_models/wall_model.lo +viscous_regularization/viscous_regularization.lo : viscous_regularization/viscous_regularization.f90 common/time_state.lo sem/dofmap.lo sem/coef.lo field/field.lo config/num_types.lo +viscous_regularization/viscous_regularization_fctry.lo : viscous_regularization/viscous_regularization_fctry.f90 common/utils.lo viscous_regularization/artificial_viscosity.lo viscous_regularization/viscous_regularization.lo +viscous_regularization/artificial_viscosity.lo : viscous_regularization/artificial_viscosity.f90 common/utils.lo sem/dofmap.lo math/field_math.lo sem/coef.lo common/json_utils.lo field/field.lo registries/registry.lo viscous_regularization/viscous_regularization.lo config/num_types.lo +viscous_regularization/avm/avm.lo : viscous_regularization/avm/avm.f90 common/time_state.lo sem/dofmap.lo sem/coef.lo registries/registry.lo field/field.lo case.lo config/num_types.lo +viscous_regularization/avm/avm_fctry.lo : viscous_regularization/avm/avm_fctry.f90 common/utils.lo viscous_regularization/avm/entropy_viscosity.lo case.lo viscous_regularization/avm/avm.lo +viscous_regularization/avm/entropy_viscosity.lo : viscous_regularization/avm/entropy_viscosity.f90 fluid/fluid_scheme_compressible.lo fluid/bcknd/device/compressible_ops_device.lo fluid/bcknd/cpu/compressible_ops_cpu.lo viscous_regularization/avm/bcknd/device/entropy_viscosity_device.lo viscous_regularization/avm/bcknd/cpu/entropy_viscosity_cpu.lo math/bcknd/device/device_math.lo device/device.lo config/neko_config.lo gs/gs_ops.lo sem/space.lo mesh/mesh.lo registries/scratch_registry.lo math/math.lo math/operators.lo time_schemes/bdf_time_scheme.lo common/time_state.lo gs/gather_scatter.lo sem/dofmap.lo sem/coef.lo field/field_series.lo math/field_math.lo field/field.lo common/json_utils.lo viscous_regularization/avm/avm.lo registries/registry.lo common/utils.lo case.lo config/num_types.lo +viscous_regularization/avm/bcknd/cpu/entropy_viscosity_cpu.lo : viscous_regularization/avm/bcknd/cpu/entropy_viscosity_cpu.f90 config/num_types.lo +viscous_regularization/avm/bcknd/device/entropy_viscosity_device.lo : viscous_regularization/avm/bcknd/device/entropy_viscosity_device.F90 common/utils.lo config/num_types.lo neko.lo : neko.f90 io/hdf5_session.lo ale/ale_manager.lo common/user_access_singleton.lo source_terms/source_term.lo common/time_step_controller.lo common/time_based_controller.lo simulation_components/lambda2.lo simulation_components/force_torque.lo simulation_components/weak_gradient_simcomp.lo simulation_components/gradient_simcomp.lo simulation_components/curl_simcomp.lo simulation_components/divergence_simcomp.lo simulation_components/derivative_simcomp.lo simulation_components/field_writer.lo les/les_model.lo common/json_utils.lo common/runtime_statistics.lo bc/field_neumann.lo bc/field_dirichlet_vector.lo bc/field_dirichlet.lo mesh/point_zone_registry.lo mesh/point_zones/sphere_point_zone.lo mesh/point_zones/box_point_zone.lo mesh/point_zone.lo sem/point_interpolator.lo common/time_interpolator.lo io/data_streamer.lo simulation_components/simcomp_executor.lo registries/scratch_registry.lo registries/registry.lo qoi/drag_torque.lo common/system.lo common/profiler.lo simulation_components/spectral_error.lo simulation_components/probes.lo common/boundary_data.lo simulation_components/boundary_flux.lo simulation_components/boundary_operation.lo simulation_components/simulation_component.lo math/tensor.lo math/matrix.lo math/vector_list.lo math/vector.lo source_terms/user_source_term.lo field/field_list.lo fluid/fluid_stats.lo sem/cpr.lo sem/map_2d.lo sem/map_1d.lo math/bcknd/device/device_math.lo device/device.lo common/jobctrl.lo common/time_state.lo common/signal.lo common/user_intf.lo common/projection.lo math/mathops.lo math/operators.lo simulation.lo io/output.lo io/output_controller.lo case.lo config/neko_config.lo comm/parmetis.lo math/ax.lo bc/dirichlet.lo bc/bc_list.lo bc/zero_dirichlet.lo bc/vector_bc_projector.lo bc/scalar_bc_projector.lo bc/mixed_bc.lo bc/bc.lo sem/coef.lo krylov/krylov.lo gs/gather_scatter.lo comm/mpi_types.lo math/field_math.lo field/field.lo io/file.lo global_interpolation/global_interpolation.lo math/mxm_wrapper.lo io/import_field_utils.lo io/format/map.lo field/mesh_field.lo mesh/point.lo mesh/mesh.lo adt/tuple.lo adt/stack.lo adt/uset.lo adt/htable.lo sem/space.lo sem/dofmap.lo sem/speclib.lo math/math.lo common/mask.lo common/log.lo common/utils.lo comm/comm.lo config/num_types.lo diff --git a/src/.depends_device b/src/.depends_device index 82fc3c3a3a16..63bb3f76b800 100644 --- a/src/.depends_device +++ b/src/.depends_device @@ -41,7 +41,7 @@ fluid/bcknd/device/hip/compressible_res.lo : fluid/bcknd/device/hip/compressible fluid/bcknd/device/hip/compressible_ops_compute_max_wave_speed.lo : fluid/bcknd/device/hip/compressible_ops_compute_max_wave_speed.hip fluid/bcknd/device/hip/compressible_ops_kernel.h fluid/bcknd/device/hip/compressible_ops_compute_entropy.lo : fluid/bcknd/device/hip/compressible_ops_compute_entropy.hip fluid/bcknd/device/hip/compressible_ops_kernel.h fluid/bcknd/device/hip/compressible_ops_update.lo : fluid/bcknd/device/hip/compressible_ops_update.hip fluid/bcknd/device/hip/compressible_ops_kernel.h -fluid/regularization/bcknd/device/hip/entropy_viscosity.lo : fluid/regularization/bcknd/device/hip/entropy_viscosity.hip fluid/regularization/bcknd/device/hip/entropy_viscosity_kernel.h +viscous_regularization/avm/bcknd/device/hip/entropy_viscosity.lo : viscous_regularization/avm/bcknd/device/hip/entropy_viscosity.hip viscous_regularization/avm/bcknd/device/hip/entropy_viscosity_kernel.h scalar/bcknd/device/hip/scalar_residual.lo : scalar/bcknd/device/hip/scalar_residual.hip scalar/bcknd/device/hip/scalar_residual_update_kernel.h sem/bcknd/device/hip/coef.lo : sem/bcknd/device/hip/coef.hip sem/bcknd/device/hip/coef_kernel.h lpt/bcknd/device/hip/lpt_periodic_bc.lo : lpt/bcknd/device/hip/lpt_periodic_bc.hip lpt/bcknd/device/hip/lpt_periodic_bc_kernel.h @@ -98,7 +98,7 @@ fluid/bcknd/device/cuda/compressible_res.lo : fluid/bcknd/device/cuda/compressib fluid/bcknd/device/cuda/compressible_ops_compute_max_wave_speed.lo : fluid/bcknd/device/cuda/compressible_ops_compute_max_wave_speed.cu fluid/bcknd/device/cuda/compressible_ops_kernel.h fluid/bcknd/device/cuda/compressible_ops_compute_entropy.lo : fluid/bcknd/device/cuda/compressible_ops_compute_entropy.cu fluid/bcknd/device/cuda/compressible_ops_kernel.h fluid/bcknd/device/cuda/compressible_ops_update.lo : fluid/bcknd/device/cuda/compressible_ops_update.cu fluid/bcknd/device/cuda/compressible_ops_kernel.h -fluid/regularization/bcknd/device/cuda/entropy_viscosity.lo : fluid/regularization/bcknd/device/cuda/entropy_viscosity.cu fluid/regularization/bcknd/device/cuda/entropy_viscosity_kernel.h +viscous_regularization/avm/bcknd/device/cuda/entropy_viscosity.lo : viscous_regularization/avm/bcknd/device/cuda/entropy_viscosity.cu viscous_regularization/avm/bcknd/device/cuda/entropy_viscosity_kernel.h fluid/stress_formulation/bcknd/device/cuda/pnpn_stress_res.lo : fluid/stress_formulation/bcknd/device/cuda/pnpn_stress_res.cu fluid/stress_formulation/bcknd/device/cuda/prs_stress_res_kernel.h common/bcknd/device/cuda/gradient_jump_penalty.lo : common/bcknd/device/cuda/gradient_jump_penalty.cu common/bcknd/device/cuda/gradient_jump_penalty_kernel.h scalar/bcknd/device/cuda/scalar_residual.lo : scalar/bcknd/device/cuda/scalar_residual.cu scalar/bcknd/device/cuda/scalar_residual_update_kernel.h @@ -150,7 +150,7 @@ fluid/bcknd/device/opencl/compressible_res.lo : fluid/bcknd/device/opencl/compre fluid/bcknd/device/opencl/compressible_ops_compute_max_wave_speed.lo : fluid/bcknd/device/opencl/compressible_ops_compute_max_wave_speed.c fluid/bcknd/device/opencl/compressible_ops_compute_max_wave_speed_kernel.cl.h fluid/bcknd/device/opencl/compressible_ops_compute_entropy.lo : fluid/bcknd/device/opencl/compressible_ops_compute_entropy.c fluid/bcknd/device/opencl/compressible_ops_compute_entropy_kernel.cl.h fluid/bcknd/device/opencl/compressible_ops_update.lo : fluid/bcknd/device/opencl/compressible_ops_update.c fluid/bcknd/device/opencl/compressible_ops_update_kernel.cl.h -fluid/regularization/bcknd/device/opencl/entropy_viscosity.lo : fluid/regularization/bcknd/device/opencl/entropy_viscosity.c fluid/regularization/bcknd/device/opencl/entropy_viscosity_kernel.cl.h +viscous_regularization/avm/bcknd/device/opencl/entropy_viscosity.lo : viscous_regularization/avm/bcknd/device/opencl/entropy_viscosity.c viscous_regularization/avm/bcknd/device/opencl/entropy_viscosity_kernel.cl.h scalar/bcknd/device/opencl/scalar_residual.lo : scalar/bcknd/device/opencl/scalar_residual.c scalar/bcknd/device/opencl/scalar_residual_kernel.cl.h sem/bcknd/device/opencl/coef.lo : sem/bcknd/device/opencl/coef.c sem/bcknd/device/opencl/coef_kernel.cl.h sem/bcknd/device/opencl/local_interpolation.lo : sem/bcknd/device/opencl/local_interpolation.c sem/bcknd/device/opencl/local_interpolation_kernel.cl.h diff --git a/src/Makefile.am b/src/Makefile.am index 9ebbb5600afe..15e7b9b4fd04 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -288,11 +288,6 @@ neko_fortran_SOURCES = \ fluid/fluid_scheme_base.f90\ fluid/fluid_scheme_incompressible.f90\ fluid/fluid_scheme_compressible.f90\ - fluid/regularization.f90\ - fluid/regularization_fctry.f90\ - fluid/regularization/entropy_viscosity.f90\ - fluid/regularization/bcknd/cpu/entropy_viscosity_cpu.f90\ - fluid/regularization/bcknd/device/entropy_viscosity_device.F90\ fluid/fluid_scheme_compressible_ns.f90\ fluid/fluid_aux.f90\ fluid/fluid_pnpn.f90\ @@ -404,6 +399,7 @@ neko_fortran_SOURCES = \ simulation_components/derivative_simcomp.f90\ simulation_components/les_simcomp.f90\ simulation_components/simulation_component_fctry.f90\ + simulation_components/avm_simcomp.f90\ source_terms/source_term.f90\ source_terms/user_source_term.f90\ source_terms/translation_source_term.f90\ @@ -488,6 +484,14 @@ neko_fortran_SOURCES = \ wall_models/bcknd/cpu/richardson_cpu.f90\ wall_models/bcknd/device/richardson_device.F90\ wall_models/wall_model_fctry.f90\ + viscous_regularization/viscous_regularization.f90\ + viscous_regularization/viscous_regularization_fctry.f90\ + viscous_regularization/artificial_viscosity.f90\ + viscous_regularization/avm/avm.f90\ + viscous_regularization/avm/avm_fctry.f90\ + viscous_regularization/avm/entropy_viscosity.f90\ + viscous_regularization/avm/bcknd/cpu/entropy_viscosity_cpu.f90\ + viscous_regularization/avm/bcknd/device/entropy_viscosity_device.F90\ neko.f90 neko_c_SOURCES = \ @@ -559,7 +563,6 @@ libneko_la_SOURCES += \ fluid/bcknd/device/hip/compressible_ops_compute_max_wave_speed.hip\ fluid/bcknd/device/hip/compressible_ops_compute_entropy.hip\ fluid/bcknd/device/hip/compressible_ops_update.hip\ - fluid/regularization/bcknd/device/hip/entropy_viscosity.hip\ fluid/stress_formulation/bcknd/device/hip/pnpn_stress_res.hip\ scalar/bcknd/device/hip/scalar_residual.hip\ sem/bcknd/device/hip/coef.hip\ @@ -580,6 +583,7 @@ libneko_la_SOURCES += \ wall_models/bcknd/device/hip/rough_log_law.hip\ wall_models/bcknd/device/hip/most.hip\ wall_models/bcknd/device/hip/richardson.hip\ + viscous_regularization/avm/bcknd/device/hip/entropy_viscosity.hip\ multigrid/bcknd/device/hip/amg_cheby.hip\ krylov/bcknd/device/hip/cheby_aux.hip endif @@ -629,7 +633,6 @@ libneko_la_SOURCES += \ fluid/bcknd/device/cuda/compressible_ops_compute_max_wave_speed.cu\ fluid/bcknd/device/cuda/compressible_ops_compute_entropy.cu\ fluid/bcknd/device/cuda/compressible_ops_update.cu\ - fluid/regularization/bcknd/device/cuda/entropy_viscosity.cu\ fluid/stress_formulation/bcknd/device/cuda/pnpn_stress_res.cu\ scalar/bcknd/device/cuda/scalar_residual.cu\ sem/bcknd/device/cuda/coef.cu\ @@ -651,6 +654,7 @@ libneko_la_SOURCES += \ wall_models/bcknd/device/cuda/rough_log_law.cu\ wall_models/bcknd/device/cuda/most.cu\ wall_models/bcknd/device/cuda/richardson.cu\ + viscous_regularization/avm/bcknd/device/cuda/entropy_viscosity.cu\ multigrid/bcknd/device/cuda/amg_cheby.cu\ krylov/bcknd/device/cuda/cheby_aux.cu @@ -736,7 +740,7 @@ libneko_la_SOURCES += \ fluid/bcknd/device/opencl/compressible_ops_compute_max_wave_speed.c\ fluid/bcknd/device/opencl/compressible_ops_compute_entropy.c\ fluid/bcknd/device/opencl/compressible_ops_update.c\ - fluid/regularization/bcknd/device/opencl/entropy_viscosity.c\ + viscous_regularization/avm/bcknd/device/opencl/entropy_viscosity.c\ scalar/bcknd/device/opencl/scalar_residual.c\ sem/bcknd/device/opencl/coef.c\ sem/bcknd/device/opencl/local_interpolation.c\ @@ -979,7 +983,7 @@ CLEANFILES += \ krylov/bcknd/device/opencl/*.cl.h\ fluid/bcknd/device/opencl/*.cl.h\ fluid/stress_formulation/bcknd/device/opencl/*.cl.h\ - fluid/regularization/bcknd/device/opencl/*.cl.h\ + viscous_regularization/avm/bcknd/device/opencl/*.cl.h\ scalar/bcknd/device/opencl/*.cl.h\ sem/bcknd/device/opencl/*.cl.h\ source_terms/bcknd/device/opencl/*.cl.h\ @@ -1042,7 +1046,6 @@ EXTRA_DIST = \ fluid/bcknd/device/cuda/vel_res_update_kernel.h\ fluid/bcknd/device/cuda/compressible_res_kernel.h\ fluid/bcknd/device/cuda/compressible_ops_kernel.h\ - fluid/regularization/bcknd/device/cuda/entropy_viscosity_kernel.h\ fluid/stress_formulation/bcknd/device/cuda/prs_stress_res_kernel.h\ scalar/bcknd/device/cuda/scalar_residual_update_kernel.h\ sem/bcknd/device/cuda/coef_kernel.h\ @@ -1061,6 +1064,7 @@ EXTRA_DIST = \ ale/bcknd/device/hip/ale_kinematics_kernel.h\ wall_models/bcknd/device/cuda/most_kernel.h\ wall_models/bcknd/device/cuda/richardson_kernel.h\ + viscous_regularization/avm/bcknd/device/cuda/entropy_viscosity_kernel.h\ bc/bcknd/device/hip/dirichlet_kernel.h\ bc/bcknd/device/hip/inflow_kernel.h\ bc/bcknd/device/hip/zero_dirichlet_kernel.h\ @@ -1102,7 +1106,6 @@ EXTRA_DIST = \ fluid/bcknd/device/hip/vel_res_update_kernel.h\ fluid/bcknd/device/hip/compressible_res_kernel.h\ fluid/bcknd/device/hip/compressible_ops_kernel.h\ - fluid/regularization/bcknd/device/hip/entropy_viscosity_kernel.h\ fluid/stress_formulation/bcknd/device/hip/prs_stress_res_kernel.h\ scalar/bcknd/device/hip/scalar_residual_update_kernel.h\ krylov/bcknd/device/cuda/pipecg_kernel.h\ @@ -1128,6 +1131,7 @@ EXTRA_DIST = \ wall_models/bcknd/device/hip/rough_log_law_kernel.h\ wall_models/bcknd/device/hip/most_kernel.h\ wall_models/bcknd/device/hip/richardson_kernel.h\ + viscous_regularization/avm/bcknd/device/hip/entropy_viscosity_kernel.h\ bc/bcknd/device/opencl/dirichlet_kernel.cl\ bc/bcknd/device/opencl/inflow_kernel.cl\ bc/bcknd/device/opencl/zero_dirichlet_kernel.cl\ @@ -1163,7 +1167,6 @@ EXTRA_DIST = \ fluid/bcknd/device/opencl/compressible_ops_compute_max_wave_speed_kernel.cl\ fluid/bcknd/device/opencl/compressible_ops_compute_entropy_kernel.cl\ fluid/bcknd/device/opencl/compressible_ops_update_kernel.cl\ - fluid/regularization/bcknd/device/opencl/entropy_viscosity_kernel.cl\ scalar/bcknd/device/opencl/scalar_residual_kernel.cl\ gs/bcknd/device/opencl/gs_kernels.cl\ krylov/bcknd/device/opencl/jacobi_kernel.cl\ @@ -1176,6 +1179,7 @@ EXTRA_DIST = \ filter/bcknd/device/cuda/mapping_kernels.h\ filter/bcknd/device/hip/mapping_kernels.h\ filter/bcknd/device/opencl/mapping_kernels.cl\ + viscous_regularization/avm/bcknd/device/opencl/entropy_viscosity_kernel.cl\ source_terms/bcknd/device/cuda/gradient_jump_penalty_kernel.h\ source_terms/bcknd/device/hip/gradient_jump_penalty_kernel.h\ math/bcknd/device/metal/math_kernel.metal\ diff --git a/src/fluid/bcknd/cpu/compressible_res_cpu.f90 b/src/fluid/bcknd/cpu/compressible_res_cpu.f90 index 032e7fd9ba67..c3b1355dbd22 100644 --- a/src/fluid/bcknd/cpu/compressible_res_cpu.f90 +++ b/src/fluid/bcknd/cpu/compressible_res_cpu.f90 @@ -78,7 +78,6 @@ module compressible_res_cpu !> @param Ax Matrix-vector product operator !> @param coef Coefficients for spatial discretization !> @param gs Gather-scatter operator for parallel communication - !> @param h Mesh size field !> @param artificial_visc Artificial viscosity field (entropy viscosity, min with low-order) !> @param mu Dynamic viscosity field (physical viscosity for momentum) !> @param kappa Thermal conductivity field (physical viscosity for energy) @@ -88,10 +87,10 @@ module compressible_res_cpu !> @param dt Time step size subroutine advance_primitive_variables_cpu(rho_field, m_x, m_y, m_z, & E, p, u, v, w, Ax, & - Ax_stress, coef, gs, h, artificial_visc, mu, kappa, bcs_vel, time, & + Ax_stress, coef, gs, artificial_visc, mu, kappa, bcs_vel, time, & rk_scheme, dt) type(field_t), intent(inout) :: rho_field, m_x, m_y, m_z, E - type(field_t), intent(in) :: p, u, v, w, h, artificial_visc, mu, kappa + type(field_t), intent(in) :: p, u, v, w, artificial_visc, mu, kappa class(Ax_t), intent(inout) :: Ax, Ax_stress type(coef_t), intent(inout) :: coef type(gs_t), intent(inout) :: gs @@ -251,7 +250,7 @@ subroutine advance_primitive_variables_cpu(rho_field, m_x, m_y, m_z, & k_E%items(i)%ptr, & temp_rho, temp_m_x, temp_m_y, temp_m_z, temp_E, & temp_p, temp_u, temp_v, temp_w, Ax, & - Ax_stress, coef, gs, h, artificial_visc, mu, kappa) + Ax_stress, coef, gs, artificial_visc, mu, kappa) end do ! Update the solution. Single parallel region covers all s stages. @@ -310,17 +309,16 @@ end subroutine advance_primitive_variables_cpu !> @param Ax Matrix-vector product operator for Laplacian terms !> @param coef Spatial discretization coefficients !> @param gs Gather-scatter operator for parallel communication - !> @param h Mesh size field !> @param artificial_visc Artificial viscosity field !> @param mu Dynamic viscosity field for physical momentum diffusion !> @param kappa Thermal conductivity field for physical energy diffusion subroutine evaluate_rhs_cpu(rhs_rho_field, rhs_m_x, rhs_m_y, rhs_m_z, rhs_E, & rho_field, m_x, m_y, m_z, E, p, u, v, w, Ax, & - Ax_stress, coef, gs, h, artificial_visc, mu, kappa) + Ax_stress, coef, gs, artificial_visc, mu, kappa) type(field_t), intent(inout) :: rhs_rho_field, & rhs_m_x, rhs_m_y, rhs_m_z, rhs_E type(field_t), intent(inout) :: rho_field, m_x, m_y, m_z, E - type(field_t), intent(in) :: p, u, v, w, h, artificial_visc, mu, kappa + type(field_t), intent(in) :: p, u, v, w, artificial_visc, mu, kappa class(Ax_t), intent(inout) :: Ax, Ax_stress type(coef_t), intent(inout) :: coef type(gs_t), intent(inout) :: gs diff --git a/src/fluid/bcknd/device/compressible_res_device.F90 b/src/fluid/bcknd/device/compressible_res_device.F90 index 0c6b50bd61ff..cbb134af98ea 100644 --- a/src/fluid/bcknd/device/compressible_res_device.F90 +++ b/src/fluid/bcknd/device/compressible_res_device.F90 @@ -1,4 +1,4 @@ -! Copyright (c) 2025, The Neko Authors +! Copyright (c) 2025-2026, The Neko Authors ! All rights reserved. ! ! Redistribution and use in source and binary forms, with or without @@ -415,10 +415,10 @@ end subroutine compressible_res_part_rk_sum_metal contains subroutine advance_primitive_variables_device(rho_field, & m_x, m_y, m_z, E, p, u, v, w, Ax, & - Ax_stress, coef, gs, h, artificial_visc, mu, kappa, bcs_vel, time, & + Ax_stress, coef, gs, artificial_visc, mu, kappa, bcs_vel, time, & rk_scheme, dt) type(field_t), intent(inout) :: rho_field, m_x, m_y, m_z, E - type(field_t), intent(in) :: p, u, v, w, h, artificial_visc, mu, kappa + type(field_t), intent(in) :: p, u, v, w, artificial_visc, mu, kappa class(Ax_t), intent(inout) :: Ax, Ax_stress type(coef_t), intent(inout) :: coef type(gs_t), intent(inout) :: gs @@ -563,7 +563,7 @@ subroutine advance_primitive_variables_device(rho_field, & k_m_y%items(i)%ptr, k_m_z%items(i)%ptr, k_E%items(i)%ptr, & temp_rho, temp_m_x, temp_m_y, temp_m_z, temp_E, & temp_p, temp_u, temp_v, temp_w, Ax, & - Ax_stress, coef, gs, h, artificial_visc, mu, kappa) + Ax_stress, coef, gs, artificial_visc, mu, kappa) end do ! Update the solution @@ -605,11 +605,11 @@ end subroutine advance_primitive_variables_device subroutine evaluate_rhs_device(rhs_rho_field, rhs_m_x, rhs_m_y, & rhs_m_z, rhs_E, rho_field, & m_x, m_y, m_z, E, p, u, v, w, Ax, & - Ax_stress, coef, gs, h, artificial_visc, mu, kappa) + Ax_stress, coef, gs, artificial_visc, mu, kappa) type(field_t), intent(inout) :: rhs_rho_field, rhs_m_x, rhs_m_y, & rhs_m_z, rhs_E type(field_t), intent(inout) :: rho_field, m_x, m_y, m_z, E - type(field_t), intent(in) :: p, u, v, w, h, artificial_visc, mu, kappa + type(field_t), intent(in) :: p, u, v, w, artificial_visc, mu, kappa class(Ax_t), intent(inout) :: Ax, Ax_stress type(coef_t), intent(inout) :: coef type(gs_t), intent(inout) :: gs diff --git a/src/fluid/compressible_res.f90 b/src/fluid/compressible_res.f90 index 91f11f53c534..0ac8bdcccca7 100644 --- a/src/fluid/compressible_res.f90 +++ b/src/fluid/compressible_res.f90 @@ -1,4 +1,4 @@ -! Copyright (c) 2025, The Neko Authors +! Copyright (c) 2025-2026, The Neko Authors ! All rights reserved. ! ! Redistribution and use in source and binary forms, with or without @@ -52,7 +52,7 @@ module compressible_residual !> Abstract interface to evaluate rhs abstract interface subroutine compressible_rhs(rho_field, m_x, m_y, m_z, E, p, u, v, w, Ax, & - Ax_stress, coef, gs, h, artificial_visc, mu, kappa, bcs_vel, time, & + Ax_stress, coef, gs, artificial_visc, mu, kappa, bcs_vel, time, & rk_scheme, dt) import field_t import Ax_t @@ -63,7 +63,7 @@ subroutine compressible_rhs(rho_field, m_x, m_y, m_z, E, p, u, v, w, Ax, & import bc_list_t import time_state_t type(field_t), intent(inout) :: rho_field, m_x, m_y, m_z, E - type(field_t), intent(in) :: p, u, v, w, h, artificial_visc, mu, kappa + type(field_t), intent(in) :: p, u, v, w, artificial_visc, mu, kappa class(Ax_t), intent(inout) :: Ax, Ax_stress type(coef_t), intent(inout) :: coef type(gs_t), intent(inout) :: gs diff --git a/src/fluid/fluid_scheme_compressible.f90 b/src/fluid/fluid_scheme_compressible.f90 index 046c230a9b28..2de1f210090f 100644 --- a/src/fluid/fluid_scheme_compressible.f90 +++ b/src/fluid/fluid_scheme_compressible.f90 @@ -1,4 +1,4 @@ -! Copyright (c) 2025, The Neko Authors +! Copyright (c) 2025-2026, The Neko Authors ! All rights reserved. ! ! Redistribution and use in source and binary forms, with or without @@ -51,11 +51,9 @@ module fluid_scheme_compressible use operators, only : cfl_compressible use device, only : device_memcpy, HOST_TO_DEVICE use compressible_ops_cpu, only : & - compressible_ops_cpu_compute_max_wave_speed, & - compressible_ops_cpu_compute_entropy + compressible_ops_cpu_compute_max_wave_speed use compressible_ops_device, only : & - compressible_ops_device_compute_max_wave_speed, & - compressible_ops_device_compute_entropy + compressible_ops_device_compute_max_wave_speed use neko_config, only : NEKO_BCKND_DEVICE use time_state, only : time_state_t use logger, only : neko_log, LOG_SIZE @@ -74,7 +72,6 @@ module fluid_scheme_compressible type(field_t), pointer :: E => null() !< Total energy type(field_t), pointer :: temperature => null() !< Temperature field type(field_t), pointer :: max_wave_speed => null() !< Maximum wave speed field - type(field_t), pointer :: S => null() !< Entropy field type(field_t), pointer :: artificial_visc => null() !< Artificial viscosity field (without physical) type(field_t), pointer :: kappa => null() !< Thermal conductivity @@ -102,9 +99,6 @@ module fluid_scheme_compressible !> Update variable material properties procedure, pass(this) :: update_material_properties => & fluid_scheme_compressible_update_material_properties - !> Compute entropy field - procedure, pass(this) :: compute_entropy => & - fluid_scheme_compressible_compute_entropy !> Compute maximum wave speed procedure, pass(this) :: compute_max_wave_speed => & fluid_scheme_compressible_compute_max_wave_speed @@ -187,11 +181,6 @@ subroutine fluid_scheme_compressible_init(this, msh, lx, params, scheme, user) this%max_wave_speed => neko_registry%get_field("max_wave_speed") call this%max_wave_speed%init(this%dm_Xh, "max_wave_speed") - ! Assign entropy field - call neko_registry%add_field(this%dm_Xh, "S") - this%S => neko_registry%get_field("S") - call this%S%init(this%dm_Xh, "S") - ! Assign artificial viscosity field (without physical viscosity) call neko_registry%add_field(this%dm_Xh, "artificial_visc") this%artificial_visc => neko_registry%get_field("artificial_visc") @@ -291,10 +280,6 @@ subroutine fluid_scheme_compressible_free(this) call this%max_wave_speed%free() end if - if (associated(this%S)) then - call this%S%free() - end if - if (associated(this%f_x)) then call this%f_x%free() deallocate(this%f_x) @@ -320,7 +305,6 @@ subroutine fluid_scheme_compressible_free(this) nullify(this%E) nullify(this%temperature) nullify(this%max_wave_speed) - nullify(this%S) nullify(this%u) nullify(this%v) @@ -458,25 +442,6 @@ subroutine fluid_scheme_compressible_update_material_properties(this, time) end subroutine fluid_scheme_compressible_update_material_properties - !> Compute entropy field S = 1/(gamma-1) * rho * (log(p) - gamma * log(rho)) - !> @param this The compressible fluid scheme object - subroutine fluid_scheme_compressible_compute_entropy(this) - class(fluid_scheme_compressible_t), intent(inout) :: this - integer :: n - - n = this%S%dof%size() - - !> TODO: Add support for SX - if (NEKO_BCKND_DEVICE .eq. 1) then - call compressible_ops_device_compute_entropy(this%S, this%p, this%rho, & - this%gamma, n) - else - call compressible_ops_cpu_compute_entropy(this%S%x, this%p%x, & - this%rho%x, this%gamma, n) - end if - - end subroutine fluid_scheme_compressible_compute_entropy - !> Compute maximum wave speed for compressible flows !> @param this The compressible fluid scheme object subroutine fluid_scheme_compressible_compute_max_wave_speed(this) @@ -509,7 +474,6 @@ subroutine fluid_scheme_compressible_log_solver_info(this, params, scheme, lx) integer, intent(in) :: lx character(len=LOG_SIZE) :: log_buf logical :: logical_val - real(kind=rp) :: real_val integer :: integer_val call neko_log%section('Fluid') @@ -538,12 +502,6 @@ subroutine fluid_scheme_compressible_log_solver_info(this, params, scheme, lx) write(log_buf, '(A,ES13.6)') 'gamma :', this%gamma call neko_log%message(log_buf) - ! Compressible-specific parameters - call json_get_or_default(params, 'case.numerics.c_avisc_low', real_val, & - 0.5_rp) - write(log_buf, '(A,ES13.6)') 'c_avisc_low:', real_val - call neko_log%message(log_buf) - call json_get_or_default(params, 'case.numerics.time_order', integer_val, 4) write(log_buf, '(A, I0)') 'RK order : ', integer_val call neko_log%message(log_buf) diff --git a/src/fluid/fluid_scheme_compressible_ns.f90 b/src/fluid/fluid_scheme_compressible_ns.f90 index 77185f8e5110..cf98d02edd94 100644 --- a/src/fluid/fluid_scheme_compressible_ns.f90 +++ b/src/fluid/fluid_scheme_compressible_ns.f90 @@ -38,12 +38,10 @@ module fluid_scheme_compressible_ns use field_series, only : field_series_t use bdf_time_scheme, only : bdf_time_scheme_t use time_scheme_controller, only : time_scheme_controller_t - use math, only : col2 - use device_math, only : device_col2 use field, only : field_t use fluid_scheme_compressible, only : fluid_scheme_compressible_t use scratch_registry, only : neko_scratch_registry - use gs_ops, only : GS_OP_ADD, GS_OP_MIN, GS_OP_MAX + use gs_ops, only : GS_OP_MIN, GS_OP_MAX use gather_scatter, only : gs_t use num_types, only : rp use mesh, only : mesh_t @@ -72,7 +70,8 @@ module fluid_scheme_compressible_ns compressible_ops_device_update_temperature use neko_config, only : NEKO_BCKND_DEVICE use mpi_f08, only : MPI_Allreduce, MPI_INTEGER, MPI_MAX - use regularization, only : regularization_t, regularization_factory + use viscous_regularization, only : viscous_regularization_t, & + viscous_regularization_factory use vector_bc_projector, only : vector_bc_projector_t, & coupled_vector_bc_projector_t implicit none @@ -82,15 +81,14 @@ module fluid_scheme_compressible_ns :: fluid_scheme_compressible_ns_t type(field_t) :: rho_res, m_x_res, m_y_res, m_z_res, m_E_res type(field_t) :: drho, dm_x, dm_y, dm_z, dE - type(field_t) :: h - real(kind=rp) :: c_avisc_low class(advection_t), allocatable :: adv class(ax_t), allocatable :: Ax class(ax_t), allocatable :: Ax_stress class(compressible_rhs_t), allocatable :: compressible_rhs type(runge_kutta_time_scheme_t) :: rk_scheme - class(regularization_t), allocatable :: regularization + class(viscous_regularization_t), allocatable :: viscous_regularization + logical :: if_viscous_regularization !> Boundary conditions projector for velocity constraints. type(coupled_vector_bc_projector_t):: bcs_vel_projector @@ -105,8 +103,7 @@ module fluid_scheme_compressible_ns !> Set up boundary conditions. procedure, pass(this) :: setup_bcs & => fluid_scheme_compressible_ns_setup_bcs - procedure, pass(this) :: compute_h - procedure, pass(this), private :: setup_regularization + procedure, pass(this), private :: setup_viscous_regularization end type fluid_scheme_compressible_ns_t interface @@ -194,7 +191,6 @@ subroutine fluid_scheme_compressible_ns_init(this, msh, lx, params, user, & call this%dm_y%init(dm_Xh, 'dm_y') call this%dm_z%init(dm_Xh, 'dm_z') call this%dE%init(dm_Xh, 'dE') - call this%h%init(dm_Xh, 'h') end associate @@ -232,11 +228,8 @@ subroutine fluid_scheme_compressible_ns_init(this, msh, lx, params, user, & ! local bases required by non-axis aligned mixed velocity conditions. call this%bcs_vel_projector%init(this%c_Xh) - ! Compute h - call this%compute_h() - - ! Initialize regularization - call this%setup_regularization(params) + ! Initialize viscous regularization + call this%setup_viscous_regularization(params) ! Initialize Runge-Kutta scheme call json_get_or_default(params, 'case.numerics.time_order', rk_order, 4) @@ -275,11 +268,10 @@ subroutine fluid_scheme_compressible_ns_free(this) call this%dm_y%free() call this%dm_z%free() call this%dE%free() - call this%h%free() - if (allocated(this%regularization)) then - call this%regularization%free() - deallocate(this%regularization) + if (allocated(this%viscous_regularization)) then + call this%viscous_regularization%free() + deallocate(this%viscous_regularization) end if do i = 1, this%bcs_density%size() @@ -300,7 +292,6 @@ end subroutine fluid_scheme_compressible_ns_free !> @param ext_bdf Time integration controller !> @param dt_controller Timestep size controller subroutine fluid_scheme_compressible_ns_step(this, time, dt_controller) - use entropy_viscosity, only : entropy_viscosity_t class(fluid_scheme_compressible_ns_t), target, intent(inout) :: this type(time_state_t), intent(in) :: time type(time_step_controller_t), intent(in) :: dt_controller @@ -320,16 +311,18 @@ subroutine fluid_scheme_compressible_ns_step(this, time, dt_controller) m_x=> this%m_x, m_y => this%m_y, m_z => this%m_z, & Xh => this%Xh, msh => this%msh, Ax => this%Ax, & c_Xh => this%c_Xh, dm_Xh => this%dm_Xh, gs_Xh => this%gs_Xh, & - E => this%E, rho => this%rho, mu => this%mu, & + E => this%E, rho => this%rho, & f_x => this%f_x, f_y => this%f_y, f_z => this%f_z, & drho => this%drho, dm_x => this%dm_x, dm_y => this%dm_y, & dm_z => this%dm_z, dE => this%dE, & - compressible_rhs => this%compressible_rhs, h => this%h, & + compressible_rhs => this%compressible_rhs, & t => time%t, tstep => time%tstep, dt => time%dt, & - c_avisc_low => this%c_avisc_low, rk_scheme => this%rk_scheme) + rk_scheme => this%rk_scheme) - ! Compute artificial viscosity - call this%regularization%compute(time) + !> Update artificial viscosity + if (allocated(this%viscous_regularization)) then + call this%viscous_regularization%update(this%artificial_visc) + end if ! Refresh user-specified physical viscosity/conductivity before RHS. call this%update_material_properties(time) @@ -337,7 +330,7 @@ subroutine fluid_scheme_compressible_ns_step(this, time, dt_controller) ! Execute RHS step with artificial viscosity field call compressible_rhs%step(rho, m_x, m_y, m_z, E, & p, u, v, w, this%Ax, & - this%Ax_stress, c_Xh, gs_Xh, h, this%artificial_visc, this%mu, & + this%Ax_stress, c_Xh, gs_Xh, this%artificial_visc, this%mu, & this%kappa, this%bcs_vel, time, rk_scheme, real(dt, kind=rp)) !> Apply density boundary conditions @@ -396,19 +389,6 @@ subroutine fluid_scheme_compressible_ns_step(this, time, dt_controller) !$omp end parallel do simd end if - !> Update entropy lag series BEFORE computing new entropy, - !> so that S_lag(1) holds the previous step's S (not the current). - !> This ensures BDF3 has 4 distinct time levels. - if (allocated(this%regularization)) then - select type (reg => this%regularization) - type is (entropy_viscosity_t) - call reg%update_lag() - end select - end if - - !> Compute entropy S = 1/(gamma-1) * rho * (log(p) - gamma * log(rho)) - call this%compute_entropy() - !> Update maximum wave speed for CFL computation call this%compute_max_wave_speed() @@ -546,77 +526,6 @@ subroutine fluid_scheme_compressible_ns_setup_bcs(this, user, params) call this%bcs_vel_projector%finalize(rebuild_mask = .true.) end subroutine fluid_scheme_compressible_ns_setup_bcs - !> Copied from les_model_compute_delta in les_model.f90 - !> TODO: move to a separate module - !> Compute characteristic mesh size h - !> @param this The fluid scheme object - subroutine compute_h(this) - class(fluid_scheme_compressible_ns_t), intent(inout) :: this - integer :: lx, ly, lz - - lx = this%c_Xh%Xh%lx - ly = this%c_Xh%Xh%ly - lz = this%c_Xh%Xh%lz - call compute_h_cpu(this%h%x, this%c_Xh%dof%x, this%c_Xh%dof%y, & - this%c_Xh%dof%z, lx, ly, lz, this%c_Xh%msh%nelv) - - if (NEKO_BCKND_DEVICE .eq. 1) then - call device_memcpy(this%h%x, this%h%x_d, this%h%dof%size(),& - HOST_TO_DEVICE, sync = .false.) - call this%gs_Xh%op(this%h, GS_OP_ADD) - call device_col2(this%h%x_d, this%c_Xh%mult_d, this%h%dof%size()) - else - call this%gs_Xh%op(this%h, GS_OP_ADD) - call col2(this%h%x, this%c_Xh%mult, this%h%dof%size()) - end if - - end subroutine compute_h - - subroutine compute_h_cpu(h, x, y, z, lx, ly, lz, nelv) - integer, intent(in) :: lx, ly, lz, nelv - real(kind=rp), intent(out) :: h(lx, ly, lz, nelv) - real(kind=rp), intent(in) :: x(lx, ly, lz, nelv) - real(kind=rp), intent(in) :: y(lx, ly, lz, nelv) - real(kind=rp), intent(in) :: z(lx, ly, lz, nelv) - integer :: e, i, j, k - integer :: im, ip, jm, jp, km, kp - real(kind=rp) :: di, dj, dk - - !$omp parallel do private(i, j, k, im, ip, jm, jp, km, kp, di, dj, dk) - do e = 1, nelv - do k = 1, lz - km = max(1, k - 1) - kp = min(lz, k + 1) - do j = 1, ly - jm = max(1, j - 1) - jp = min(ly, j + 1) - do i = 1, lx - im = max(1, i - 1) - ip = min(lx, i + 1) - - di = (x(ip, j, k, e) - x(im, j, k, e))**2 & - + (y(ip, j, k, e) - y(im, j, k, e))**2 & - + (z(ip, j, k, e) - z(im, j, k, e))**2 - - dj = (x(i, jp, k, e) - x(i, jm, k, e))**2 & - + (y(i, jp, k, e) - y(i, jm, k, e))**2 & - + (z(i, jp, k, e) - z(i, jm, k, e))**2 - - dk = (x(i, j, kp, e) - x(i, j, km, e))**2 & - + (y(i, j, kp, e) - y(i, j, km, e))**2 & - + (z(i, j, kp, e) - z(i, j, km, e))**2 - - di = sqrt(di) / (ip - im) - dj = sqrt(dj) / (jp - jm) - dk = sqrt(dk) / (kp - km) - h(i,j,k,e) = (di * dj * dk)**(1.0_rp / 3.0_rp) - end do - end do - end do - end do - !$omp end parallel do - end subroutine compute_h_cpu - !> Restart the simulation from saved state !! @param this The fluid scheme object !! @param dtlag Previous timestep sizes @@ -626,46 +535,38 @@ subroutine fluid_scheme_compressible_ns_restart(this, chkp) type(chkp_t), intent(inout) :: chkp end subroutine fluid_scheme_compressible_ns_restart - subroutine setup_regularization(this, params) - use entropy_viscosity, only : entropy_viscosity_t, & - entropy_viscosity_set_fields + subroutine setup_viscous_regularization(this, params) class(fluid_scheme_compressible_ns_t), target, intent(inout) :: this type(json_file), intent(inout) :: params type(json_file) :: reg_json - type(json_core) :: json_core_inst - type(json_value), pointer :: reg_params - character(len=:), allocatable :: buffer - real(kind=rp) :: c_avisc_entropy_val - character(len=:), allocatable :: regularization_type - - call json_get_or_default(params, 'case.numerics.c_avisc_low', & - this%c_avisc_low, 0.5_rp) - call json_get_or_default(params, 'case.numerics.c_avisc_entropy', & - c_avisc_entropy_val, 1.0_rp) - - call json_core_inst%initialize() - call json_core_inst%create_object(reg_params, '') - call json_core_inst%add(reg_params, 'c_avisc_entropy', c_avisc_entropy_val) - call json_core_inst%add(reg_params, 'c_avisc_low', this%c_avisc_low) - call json_core_inst%print_to_string(reg_params, buffer) - call json_core_inst%destroy(reg_params) - - call reg_json%initialize() - call reg_json%load_from_string(buffer) - - regularization_type = 'entropy_viscosity' - - call regularization_factory(this%regularization, regularization_type, & - reg_json, this%c_Xh, this%dm_Xh, this%artificial_visc) - - select type (reg => this%regularization) - type is (entropy_viscosity_t) - call entropy_viscosity_set_fields(reg, this%S, this%u, this%v, this%w, & - this%h, this%max_wave_speed, this%msh, this%Xh, this%gs_Xh) - end select + logical :: found + character(len=:), allocatable :: viscous_regularization_type + + found = .false. + if (this%params%valid_path('case.fluid.viscous_regularization')) then + call json_get(params, 'case.fluid.viscous_regularization', & + reg_json) + found = .true. + end if + + if (.not. found) then + ! No viscous_regularization specified, so we skip setup + this%if_viscous_regularization = .false. + return + else + this%if_viscous_regularization = .true. + end if + + call json_get(reg_json, 'type', viscous_regularization_type) + call viscous_regularization_factory(this%viscous_regularization, & + viscous_regularization_type, & + reg_json, this%c_Xh, this%dm_Xh) call reg_json%destroy() - end subroutine setup_regularization + if (allocated(viscous_regularization_type)) & + deallocate(viscous_regularization_type) + + end subroutine setup_viscous_regularization end module fluid_scheme_compressible_ns diff --git a/src/fluid/regularization.f90 b/src/fluid/regularization.f90 deleted file mode 100644 index 4963f58227db..000000000000 --- a/src/fluid/regularization.f90 +++ /dev/null @@ -1,119 +0,0 @@ -! Copyright (c) 2025, The Neko Authors -! All rights reserved. -! -! Redistribution and use in source and binary forms, with or without -! modification, are permitted provided that the following conditions -! are met: -! -! * Redistributions of source code must retain the above copyright -! notice, this list of conditions and the following disclaimer. -! -! * Redistributions in binary form must reproduce the above -! copyright notice, this list of conditions and the following -! disclaimer in the documentation and/or other materials provided -! with the distribution. -! -! * Neither the name of the authors nor the names of its -! contributors may be used to endorse or promote products derived -! from this software without specific prior written permission. -! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, -! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, -! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; -! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT -! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN -! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -! POSSIBILITY OF SUCH DAMAGE. -! -module regularization - use num_types, only : rp - use json_module, only : json_file - use field, only : field_t - use coefs, only : coef_t - use dofmap, only : dofmap_t - use time_state, only : time_state_t - implicit none - private - - type, abstract, public :: regularization_t - type(field_t), pointer :: reg_coeff => null() - type(coef_t), pointer :: coef => null() - type(dofmap_t), pointer :: dof => null() - contains - procedure, pass(this) :: init_base => regularization_init_base - procedure, pass(this) :: free_base => regularization_free_base - procedure(reg_init), pass(this), deferred :: init - procedure(reg_free), pass(this), deferred :: free - procedure(reg_compute), pass(this), deferred :: compute - end type regularization_t - - abstract interface - subroutine reg_init(this, json, coef, dof, reg_coeff) - import regularization_t, json_file, coef_t, dofmap_t, field_t - class(regularization_t), intent(inout) :: this - type(json_file), intent(inout) :: json - type(coef_t), intent(in), target :: coef - type(dofmap_t), intent(in), target :: dof - type(field_t), intent(in), target :: reg_coeff - end subroutine reg_init - end interface - - abstract interface - subroutine reg_free(this) - import regularization_t - class(regularization_t), intent(inout) :: this - end subroutine reg_free - end interface - - abstract interface - subroutine reg_compute(this, time) - import regularization_t, time_state_t - class(regularization_t), intent(inout) :: this - type(time_state_t), intent(in) :: time - end subroutine reg_compute - end interface - - interface - module subroutine regularization_factory(object, type_name, json, & - coef, dof, reg_coeff) - class(regularization_t), allocatable, intent(inout) :: object - character(len=*), intent(in) :: type_name - type(json_file), intent(inout) :: json - type(coef_t), intent(in), target :: coef - type(dofmap_t), intent(in), target :: dof - type(field_t), intent(in), target :: reg_coeff - end subroutine regularization_factory - end interface - - public :: regularization_factory - -contains - - subroutine regularization_init_base(this, json, coef, dof, reg_coeff) - class(regularization_t), intent(inout) :: this - type(json_file), intent(inout) :: json - type(coef_t), intent(in), target :: coef - type(dofmap_t), intent(in), target :: dof - type(field_t), intent(in), target :: reg_coeff - - this%coef => coef - this%dof => dof - this%reg_coeff => reg_coeff - - end subroutine regularization_init_base - - subroutine regularization_free_base(this) - class(regularization_t), intent(inout) :: this - - nullify(this%coef) - nullify(this%dof) - nullify(this%reg_coeff) - - end subroutine regularization_free_base - -end module regularization diff --git a/src/simulation_components/avm_simcomp.f90 b/src/simulation_components/avm_simcomp.f90 new file mode 100644 index 000000000000..7c9517f8be67 --- /dev/null +++ b/src/simulation_components/avm_simcomp.f90 @@ -0,0 +1,148 @@ +! Copyright (c) 2026, The Neko Authors +! All rights reserved. +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! +! * Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! +! * Redistributions in binary form must reproduce the above +! copyright notice, this list of conditions and the following +! disclaimer in the documentation and/or other materials provided +! with the distribution. +! +! * Neither the name of the authors nor the names of its +! contributors may be used to endorse or promote products derived +! from this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE +! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +!> Implements the `avm_simcomp_t` type. + +module avm_simcomp + use num_types, only : rp + use json_module, only : json_file + use simulation_component, only : simulation_component_t + use case, only : case_t + use time_state, only : time_state_t + use artificial_viscosity_model, only : avm_t, avm_factory + use json_utils, only : json_get, json_get_or_default + use field_writer, only : field_writer_t + use utils, only : neko_error, NEKO_VARNAME_LEN + implicit none + private + + !> A simulation component that drives the computation of the artificial + !! viscosity method. + type, public, extends(simulation_component_t) :: avm_simcomp_t + !> The artificial viscosity model. + class(avm_t), allocatable :: avm + !> Output writer. + type(field_writer_t) :: writer + contains + !> Constructor from json, wrapping the actual constructor. + procedure, pass(this) :: init => avm_simcomp_init_from_json + !> Destructor. + procedure, pass(this) :: free => avm_simcomp_free + !> Compute the avm_simcomp field before the time step. + procedure, pass(this) :: preprocess_ => avm_simcomp_preprocess + !> Compute the avm_simcomp field. + procedure, pass(this) :: compute_ => avm_simcomp_compute + !> Compute the avm_simcomp field when restart. + procedure, pass(this) :: restart_ => avm_simcomp_restart + end type avm_simcomp_t + +contains + + !> Constructor from json. + !! @param this The artificial viscosity simulation component. + !! @param json The simulation component configuration. + !! @param case The case being initialized. + subroutine avm_simcomp_init_from_json(this, json, case) + class(avm_simcomp_t), intent(inout), target :: this + type(json_file), intent(inout) :: json + class(case_t), intent(inout), target :: case + character(len=:), allocatable :: name + character(len=:), allocatable :: model_name + character(len=NEKO_VARNAME_LEN) :: fields(1) + + call this%free() + + call json_get_or_default(json, "name", name, "artificial_viscosity_model") + call json_get(json, "model", model_name) + this%name = name + + call this%init_base(json, case) + + ! Create the AVM first so we can get the field name it uses + call avm_factory(this%avm, model_name, case, json) + + ! Get the field name from the model's reg_coeff field and add it to the + ! field_writer output list + fields(1) = this%avm%reg_coeff%name + call json%add("fields", fields) + + call this%writer%init(json, case) + + if (allocated(name)) deallocate(name) + if (allocated(model_name)) deallocate(model_name) + end subroutine avm_simcomp_init_from_json + + !> Destructor. + !! @param this The artificial viscosity simulation component. + subroutine avm_simcomp_free(this) + class(avm_simcomp_t), intent(inout) :: this + call this%free_base() + call this%writer%free() + + if (allocated(this%avm)) then + call this%avm%free() + deallocate(this%avm) + end if + end subroutine avm_simcomp_free + + !> Compute the avm_simcomp field. + !! @param this The artificial viscosity simulation component. + !! @param time The current time info + subroutine avm_simcomp_preprocess(this, time) + class(avm_simcomp_t), intent(inout) :: this + type(time_state_t), intent(in) :: time + + call this%avm%preprocess(time) + end subroutine avm_simcomp_preprocess + + !> Compute the avm_simcomp field. + !! @param this The artificial viscosity simulation component. + !! @param time The current time info + subroutine avm_simcomp_compute(this, time) + class(avm_simcomp_t), intent(inout) :: this + type(time_state_t), intent(in) :: time + + call this%avm%compute(time) + end subroutine avm_simcomp_compute + + !> Compute the avm_simcomp field when restart. + !! @param this The artificial viscosity simulation component. + !! @param time The current time info + subroutine avm_simcomp_restart(this, time) + class(avm_simcomp_t), intent(inout) :: this + type(time_state_t), intent(in) :: time + + call this%avm%restart(time) + end subroutine avm_simcomp_restart + +end module avm_simcomp diff --git a/src/simulation_components/simulation_component_fctry.f90 b/src/simulation_components/simulation_component_fctry.f90 index b79cd161f546..b56f9f7d5bbe 100644 --- a/src/simulation_components/simulation_component_fctry.f90 +++ b/src/simulation_components/simulation_component_fctry.f90 @@ -56,13 +56,14 @@ use spectral_error, only : spectral_error_t use data_streamer_simcomp, only : data_streamer_simcomp_t use field_subsampler, only : field_subsampler_t + use avm_simcomp, only : avm_simcomp_t use wall_shear_stress_simcomp, only : wall_shear_stress_t use boundary_data_writer_simcomp, only : boundary_data_writer_t use utils, only : neko_type_error, neko_type_registration_error implicit none ! List of all possible types created by the factory routine - character(len=23) :: SIMCOMPS_KNOWN_TYPES(24) = [character(len=23) :: & + character(len=30) :: SIMCOMPS_KNOWN_TYPES(25) = [character(len=30) :: & "boundary_operation", & "boundary_flux", & "lagrangian_particles", & @@ -85,6 +86,7 @@ "spectral_error", & "data_streamer", & "field_subsampler", & + "artificial_viscosity_model", & "wall_shear_stress", & "boundary_data_writer"] @@ -175,6 +177,8 @@ module subroutine simulation_component_allocator(object, type_name) allocate(data_streamer_simcomp_t::object) case ("field_subsampler") allocate(field_subsampler_t::object) + case ("artificial_viscosity_model") + allocate(avm_simcomp_t::object) case ("wall_shear_stress") allocate(wall_shear_stress_t::object) case ("boundary_data_writer") diff --git a/src/viscous_regularization/artificial_viscosity.f90 b/src/viscous_regularization/artificial_viscosity.f90 new file mode 100644 index 000000000000..578fee76355a --- /dev/null +++ b/src/viscous_regularization/artificial_viscosity.f90 @@ -0,0 +1,109 @@ +! Copyright (c) 2026, The Neko Authors +! All rights reserved. +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! +! * Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! +! * Redistributions in binary form must reproduce the above +! copyright notice, this list of conditions and the following +! disclaimer in the documentation and/or other materials provided +! with the distribution. +! +! * Neither the name of the authors nor the names of its +! contributors may be used to endorse or promote products derived +! from this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE +! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +!> Implements `artificial_viscosity_t`. +module artificial_viscosity + use num_types, only : rp + use viscous_regularization, only : viscous_regularization_t + use registry, only : neko_registry + use field, only : field_t + use json_module, only : json_file + use json_utils, only : json_get + use coefs, only : coef_t + use field_math, only : field_copy, field_add2 + use dofmap, only : dofmap_t + use utils, only : neko_error + implicit none + private + + !> The handler for the artificial viscosity. + type, public, extends(viscous_regularization_t) :: artificial_viscosity_t + contains + !> Initialize the artificial viscosity regularization. + procedure, pass(this) :: init => artificial_viscosity_init_from_json + !> Free the artificial viscosity regularization. + procedure, pass(this) :: free => artificial_viscosity_free + !> Update the effective viscosity field. + procedure, pass(this) :: update => artificial_viscosity_update + end type artificial_viscosity_t + +contains + !> Construct an artificial viscosity regularization from JSON. + !! @param this The artificial viscosity regularization. + !! @param json The regularization configuration. + !! @param coef The SEM coefficients. + !! @param dof The SEM map of degrees of freedom. + subroutine artificial_viscosity_init_from_json(this, json, coef, dof) + class(artificial_viscosity_t), intent(inout) :: this + type(json_file), intent(inout) :: json + type(coef_t), intent(in), target :: coef + type(dofmap_t), intent(in), target :: dof + + call this%free() + call this%init_base(json, coef, dof) + + call json_get(json, "reg_coeff_name", & + this%reg_coeff_name) + + call neko_registry%add_field(dof, this%reg_coeff_name) + this%reg_coeff => neko_registry%get_field(this%reg_coeff_name) + + end subroutine artificial_viscosity_init_from_json + + !> Free an artificial viscosity regularization. + !! @param this The artificial viscosity regularization. + subroutine artificial_viscosity_free(this) + class(artificial_viscosity_t), intent(inout) :: this + + nullify(this%reg_coeff) + call this%free_base() + + end subroutine artificial_viscosity_free + + !> Update the effective viscosity with artificial and physical viscosity. + !! @param this The artificial viscosity regularization. + !! @param effective_visc The effective viscosity field to update. + !! @param mu The optional physical viscosity field. + subroutine artificial_viscosity_update(this, effective_visc, mu) + class(artificial_viscosity_t), intent(inout) :: this + type(field_t), intent(inout) :: effective_visc + type(field_t), intent(in), optional :: mu + + call field_copy(effective_visc, this%reg_coeff, this%coef%dof%size()) + if (present(mu)) then + call field_add2(effective_visc, mu, this%coef%dof%size()) + end if + + end subroutine artificial_viscosity_update + +end module artificial_viscosity diff --git a/src/viscous_regularization/avm/avm.f90 b/src/viscous_regularization/avm/avm.f90 new file mode 100644 index 000000000000..fc7a1e385106 --- /dev/null +++ b/src/viscous_regularization/avm/avm.f90 @@ -0,0 +1,235 @@ +! Copyright (c) 2026, The Neko Authors +! All rights reserved. +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! +! * Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! +! * Redistributions in binary form must reproduce the above +! copyright notice, this list of conditions and the following +! disclaimer in the documentation and/or other materials provided +! with the distribution. +! +! * Neither the name of the authors nor the names of its +! contributors may be used to endorse or promote products derived +! from this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE +! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +!> Implements `avm_t`. +module artificial_viscosity_model + use num_types, only : rp + use case, only : case_t + use field, only : field_t + use json_module, only : json_file + use registry, only : neko_registry + use coefs, only : coef_t + use dofmap, only : dofmap_t + use time_state, only : time_state_t + implicit none + private + + !> Base abstract type for artificial viscosity models. + type, abstract, public :: avm_t + !> Artificial viscosity. + type(field_t), pointer :: reg_coeff => null() + !> SEM coefficients. + type(coef_t), pointer :: coef => null() + !> SEM dofmap. + type(dofmap_t), pointer :: dof => null() + !> Number of custom artificial viscosity model allocators. + integer :: les_model_registry_size = 0 + contains + !> Constructor for the avm_t (base) class. + procedure, pass(this) :: init_base => avm_init_base + !> Destructor for the avm_t (base) class. + procedure, pass(this) :: free_base => avm_free_base + !> The common constructor. + procedure(avm_init), pass(this), deferred :: init + !> Destructor. + procedure(avm_free), pass(this), deferred :: free + !> Perform artificial viscosity related computations before the time step. + procedure(avm_preprocess), pass(this), deferred :: preprocess + !> Perform artificial viscosity related computations after the time step. + procedure(avm_compute), pass(this), deferred :: compute + !> Restore model state after loading a checkpoint. + procedure, pass(this) :: restart => avm_restart + end type avm_t + + abstract interface + !> Perform artificial viscosity related computations before the time step. + !! @param this The artificial viscosity model. + !! @param time The time state. + subroutine avm_preprocess(this, time) + import avm_t, time_state_t + class(avm_t), intent(inout) :: this + type(time_state_t), intent(in) :: time + end subroutine avm_preprocess + end interface + + abstract interface + !> Perform artificial viscosity related computations after the time step. + !! @param this The artificial viscosity model. + !! @param time The time state. + subroutine avm_compute(this, time) + import avm_t, time_state_t + class(avm_t), intent(inout) :: this + type(time_state_t), intent(in) :: time + end subroutine avm_compute + end interface + + abstract interface + !> Common constructor. + !! @param this The artificial viscosity model. + !! @param case The case_t object. + !! @param json A dictionary with parameters. + subroutine avm_init(this, case, json) + import avm_t, json_file, case_t + class(avm_t), intent(inout) :: this + class(case_t), intent(inout), target :: case + type(json_file), intent(inout) :: json + end subroutine avm_init + end interface + + abstract interface + !> Destructor. + !! @param this The artificial viscosity model. + subroutine avm_free(this) + import avm_t + class(avm_t), intent(inout) :: this + end subroutine avm_free + end interface + + interface + !> Artificial viscosity model allocator. + !! @param object The object to be allocated. + !! @param type_name The name of the artificial viscosity model. + module subroutine avm_allocator(object, type_name) + class(avm_t), allocatable, intent(inout) :: object + character(len=*), intent(in) :: type_name + end subroutine avm_allocator + end interface + + ! + ! Machinery for injecting user-defined types + ! + + !> Interface for an object allocator. + !! Implemented in the user modules, should allocate the `obj` to the custom + !! user type. + abstract interface + !> Allocate a custom artificial viscosity model. + !! @param obj The object to allocate. + subroutine avm_allocate(obj) + import avm_t + class(avm_t), allocatable, intent(inout) :: obj + end subroutine avm_allocate + end interface + + interface + !> Called in user modules to add an allocator for custom types. + !! @param type_name The custom model type name. + !! @param allocator The custom model allocator. + module subroutine register_avm(type_name, allocator) + character(len=*), intent(in) :: type_name + procedure(avm_allocate), pointer, intent(in) :: allocator + end subroutine register_avm + end interface + + !> A name-allocator pair for user-defined types. A helper type to define a + !! registry of custom allocators. + type allocator_entry + !> Registered model type name. + character(len=20) :: type_name + !> Allocator associated with the model type. + procedure(avm_allocate), pointer, nopass :: allocator + end type allocator_entry + + !> Registry of artificial viscosity model allocators for user-defined types + type(allocator_entry), allocatable :: avm_registry(:) + + !> The size of the `avm_registry` + integer :: avm_registry_size = 0 + + + + interface + !> Artificial viscosity model factory. Both allocates and initializes + !! the object. + !! @param object The object to be allocated. + !! @param type_name The name of the artificial viscosity model. + !! @param case The case_t object. + !! @param json A dictionary with parameters. + module subroutine avm_factory(object, type_name, case, json) + class(avm_t), allocatable, intent(inout) :: object + character(len=*), intent(in) :: type_name + class(case_t), intent(inout), target :: case + type(json_file), intent(inout) :: json + end subroutine avm_factory + end interface + + public :: avm_factory, avm_allocator, register_avm, & + avm_allocate + + +contains + !> Constructor for the avm_t (base) class. + !! @param this The artificial viscosity model. + !! @param dof Map of degrees of freedom. + !! @param coef The SEM coefficients. + !! @param reg_coeff_name The name of the artificial viscosity field. + subroutine avm_init_base(this, dof, coef, reg_coeff_name) + class(avm_t), intent(inout) :: this + class(dofmap_t), intent(in), target :: dof + class(coef_t), intent(in), target :: coef + character(len=*), intent(in) :: reg_coeff_name + + call this%free_base() + + call neko_registry%add_field(dof, trim(reg_coeff_name), .true.) + + this%reg_coeff => neko_registry%get_field(trim(reg_coeff_name)) + this%coef => coef + this%dof => dof + + end subroutine avm_init_base + + !> Restore model state after loading a checkpoint. + !! Models with additional history should override this method. + !! @param this The artificial viscosity model. + !! @param time The restored time state. + subroutine avm_restart(this, time) + class(avm_t), intent(inout) :: this + type(time_state_t), intent(in) :: time + + call this%compute(time) + + end subroutine avm_restart + + !> Destructor for the avm_t (base) class. + !! @param this The artificial viscosity model. + subroutine avm_free_base(this) + class(avm_t), intent(inout) :: this + + nullify(this%reg_coeff) + nullify(this%coef) + nullify(this%dof) + + end subroutine avm_free_base + +end module artificial_viscosity_model diff --git a/src/viscous_regularization/avm/avm_fctry.f90 b/src/viscous_regularization/avm/avm_fctry.f90 new file mode 100644 index 000000000000..dda49ea898e7 --- /dev/null +++ b/src/viscous_regularization/avm/avm_fctry.f90 @@ -0,0 +1,130 @@ +! Copyright (c) 2026, The Neko Authors +! All rights reserved. +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! +! * Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! +! * Redistributions in binary form must reproduce the above +! copyright notice, this list of conditions and the following +! disclaimer in the documentation and/or other materials provided +! with the distribution. +! +! * Neither the name of the authors nor the names of its +! contributors may be used to endorse or promote products derived +! from this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE +! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +!> Implements the artificial viscosity model factory and allocator. +submodule (artificial_viscosity_model) artificial_viscosity_model_fctry + use case, only : case_t + use entropy_viscosity, only : entropy_viscosity_t + use utils, only : neko_type_error, neko_type_registration_error + implicit none + + !> List of all built-in artificial viscosity model types. created by the factory routine + character(len=20) :: AVM_KNOWN_TYPES(1) = [character(len=20) :: & + "entropy_viscosity"] + +contains + !> artificial viscosity model factory. + !! @param object The object to be allocated. + !! @param type_name The name of the artificial viscosity model. + !! @param case The case_t object. + !! @param json A dictionary with parameters. + module subroutine avm_factory(object, type_name, case, json) + class(avm_t), allocatable, intent(inout) :: object + character(len=*), intent(in) :: type_name + class(case_t), intent(inout), target :: case + type(json_file), intent(inout) :: json + character(len=:), allocatable :: type_string + + call avm_allocator(object, type_name) + call object%init(case, json) + end subroutine avm_factory + + !> artificial viscosity model allocator. + !! @param object The object to be allocated. + !! @param type_name The name of the artificial viscosity model. + module subroutine avm_allocator(object, type_name) + class(avm_t), allocatable, intent(inout) :: object + character(len=*), intent(in) :: type_name + integer :: i + + if (allocated(object)) then + call object%free() + deallocate(object) + end if + + select case (trim(type_name)) + case ('entropy_viscosity') + allocate(entropy_viscosity_t::object) + case default + do i = 1, avm_registry_size + if (trim(type_name) == trim(avm_registry(i)%type_name)) then + call avm_registry(i)%allocator(object) + return + end if + end do + + call neko_type_error("artificial viscosity model", & + type_name, AVM_KNOWN_TYPES) + end select + + end subroutine avm_allocator + + !> Register a custom artificial viscosity model allocator. + !! Called in custom user module inside the `module_name_register_types` + !! routine to add a custom type allocator to the registry. + !! @param type_name The name of the type to allocate. + !! @param allocator The allocator for the custom user type. + module subroutine register_avm(type_name, allocator) + character(len=*), intent(in) :: type_name + procedure(avm_allocate), pointer, intent(in) :: allocator + type(allocator_entry), allocatable :: temp(:) + integer :: i + + do i = 1, size(AVM_KNOWN_TYPES) + if (trim(type_name) .eq. trim(AVM_KNOWN_TYPES(i))) then + call neko_type_registration_error("artificial viscosity model", & + type_name, .true.) + end if + end do + + do i = 1, avm_registry_size + if (trim(type_name) .eq. trim(avm_registry(i)%type_name)) then + call neko_type_registration_error("artificial viscosity model", & + type_name, .false.) + end if + end do + + ! Expand registry + if (avm_registry_size == 0) then + allocate(avm_registry(1)) + else + allocate(temp(avm_registry_size + 1)) + temp(1:avm_registry_size) = avm_registry + call move_alloc(temp, avm_registry) + end if + + avm_registry_size = avm_registry_size + 1 + avm_registry(avm_registry_size)%type_name = type_name + avm_registry(avm_registry_size)%allocator => allocator + end subroutine register_avm + +end submodule artificial_viscosity_model_fctry diff --git a/src/fluid/regularization/bcknd/cpu/entropy_viscosity_cpu.f90 b/src/viscous_regularization/avm/bcknd/cpu/entropy_viscosity_cpu.f90 similarity index 99% rename from src/fluid/regularization/bcknd/cpu/entropy_viscosity_cpu.f90 rename to src/viscous_regularization/avm/bcknd/cpu/entropy_viscosity_cpu.f90 index 6d4dcf9884d6..faeef9172b28 100644 --- a/src/fluid/regularization/bcknd/cpu/entropy_viscosity_cpu.f90 +++ b/src/viscous_regularization/avm/bcknd/cpu/entropy_viscosity_cpu.f90 @@ -1,4 +1,4 @@ -! Copyright (c) 2025, The Neko Authors +! Copyright (c) 2025-2026, The Neko Authors ! All rights reserved. ! ! Redistribution and use in source and binary forms, with or without diff --git a/src/fluid/regularization/bcknd/device/cuda/entropy_viscosity.cu b/src/viscous_regularization/avm/bcknd/device/cuda/entropy_viscosity.cu similarity index 100% rename from src/fluid/regularization/bcknd/device/cuda/entropy_viscosity.cu rename to src/viscous_regularization/avm/bcknd/device/cuda/entropy_viscosity.cu diff --git a/src/fluid/regularization/bcknd/device/cuda/entropy_viscosity_kernel.h b/src/viscous_regularization/avm/bcknd/device/cuda/entropy_viscosity_kernel.h similarity index 100% rename from src/fluid/regularization/bcknd/device/cuda/entropy_viscosity_kernel.h rename to src/viscous_regularization/avm/bcknd/device/cuda/entropy_viscosity_kernel.h diff --git a/src/fluid/regularization/bcknd/device/entropy_viscosity_device.F90 b/src/viscous_regularization/avm/bcknd/device/entropy_viscosity_device.F90 similarity index 90% rename from src/fluid/regularization/bcknd/device/entropy_viscosity_device.F90 rename to src/viscous_regularization/avm/bcknd/device/entropy_viscosity_device.F90 index 388e122754e7..5d90dfcf05c4 100644 --- a/src/fluid/regularization/bcknd/device/entropy_viscosity_device.F90 +++ b/src/viscous_regularization/avm/bcknd/device/entropy_viscosity_device.F90 @@ -1,4 +1,4 @@ -! Copyright (c) 2025, The Neko Authors +! Copyright (c) 2025-2026, The Neko Authors ! All rights reserved. ! ! Redistribution and use in source and binary forms, with or without @@ -278,7 +278,15 @@ end subroutine metal_entropy_visc_smooth_divide contains - !> Compute entropy residual on device + !> Compute entropy residual on a device. + !! @param entropy_residual_d Device entropy residual field. + !! @param S_d Device current entropy field. + !! @param S_lag1_d Device first lagged entropy field. + !! @param S_lag2_d Device second lagged entropy field. + !! @param S_lag3_d Device third lagged entropy field. + !! @param bdf_coeffs BDF time-scheme coefficients. + !! @param dt Time-step size. + !! @param n Number of points. subroutine entropy_viscosity_compute_residual_device(entropy_residual_d, & S_d, S_lag1_d, S_lag2_d, S_lag3_d, bdf_coeffs, dt, n) type(c_ptr), intent(in) :: entropy_residual_d @@ -308,7 +316,13 @@ subroutine entropy_viscosity_compute_residual_device(entropy_residual_d, & #endif end subroutine entropy_viscosity_compute_residual_device - !> Compute viscosity from entropy residual on device + !> Compute viscosity from an entropy residual on a device. + !! @param reg_coeff_d Device regularization coefficient field. + !! @param entropy_residual_d Device entropy residual field. + !! @param h_d Device mesh-size field. + !! @param c_avisc_entropy Entropy viscosity coefficient. + !! @param n_S Entropy normalization factor. + !! @param n Number of points. subroutine entropy_viscosity_compute_viscosity_device(reg_coeff_d, & entropy_residual_d, h_d, c_avisc_entropy, n_S, n) type(c_ptr), intent(in) :: reg_coeff_d, entropy_residual_d, h_d @@ -332,7 +346,10 @@ subroutine entropy_viscosity_compute_viscosity_device(reg_coeff_d, & #endif end subroutine entropy_viscosity_compute_viscosity_device - !> Apply element-wise maximum on device + !> Apply the element-wise maximum on a device. + !! @param reg_coeff_d Device regularization coefficient field. + !! @param lx Number of points per element direction. + !! @param nelv Number of local elements. subroutine entropy_viscosity_apply_element_max_device(reg_coeff_d, lx, nelv) type(c_ptr), intent(in) :: reg_coeff_d integer, intent(in) :: lx, nelv @@ -350,7 +367,12 @@ subroutine entropy_viscosity_apply_element_max_device(reg_coeff_d, lx, nelv) #endif end subroutine entropy_viscosity_apply_element_max_device - !> Clamp regularization coefficient to low-order viscosity on device + !> Clamp the coefficient to low-order viscosity on a device. + !! @param reg_coeff_d Device regularization coefficient field. + !! @param h_d Device mesh-size field. + !! @param max_wave_speed_d Device maximum wave-speed field. + !! @param c_avisc_low Low-order viscosity coefficient. + !! @param n Number of points. subroutine entropy_viscosity_clamp_to_low_order_device(reg_coeff_d, & h_d, max_wave_speed_d, c_avisc_low, n) type(c_ptr), intent(in) :: reg_coeff_d, h_d, max_wave_speed_d @@ -374,7 +396,11 @@ subroutine entropy_viscosity_clamp_to_low_order_device(reg_coeff_d, & #endif end subroutine entropy_viscosity_clamp_to_low_order_device - !> Divide by multiplicity for smoothing on device + !> Divide by gather-scatter multiplicity on a device. + !! @param reg_coeff_d Device regularization coefficient field. + !! @param temp_field_d Device field containing summed values. + !! @param mult_field_d Device multiplicity field. + !! @param n Number of points. subroutine entropy_viscosity_smooth_divide_device(reg_coeff_d, & temp_field_d, mult_field_d, n) type(c_ptr), intent(in) :: reg_coeff_d, temp_field_d, mult_field_d diff --git a/src/fluid/regularization/bcknd/device/hip/entropy_viscosity.hip b/src/viscous_regularization/avm/bcknd/device/hip/entropy_viscosity.hip similarity index 100% rename from src/fluid/regularization/bcknd/device/hip/entropy_viscosity.hip rename to src/viscous_regularization/avm/bcknd/device/hip/entropy_viscosity.hip diff --git a/src/fluid/regularization/bcknd/device/hip/entropy_viscosity_kernel.h b/src/viscous_regularization/avm/bcknd/device/hip/entropy_viscosity_kernel.h similarity index 100% rename from src/fluid/regularization/bcknd/device/hip/entropy_viscosity_kernel.h rename to src/viscous_regularization/avm/bcknd/device/hip/entropy_viscosity_kernel.h diff --git a/src/fluid/regularization/bcknd/device/opencl/entropy_viscosity.c b/src/viscous_regularization/avm/bcknd/device/opencl/entropy_viscosity.c similarity index 100% rename from src/fluid/regularization/bcknd/device/opencl/entropy_viscosity.c rename to src/viscous_regularization/avm/bcknd/device/opencl/entropy_viscosity.c diff --git a/src/fluid/regularization/bcknd/device/opencl/entropy_viscosity_kernel.cl b/src/viscous_regularization/avm/bcknd/device/opencl/entropy_viscosity_kernel.cl similarity index 100% rename from src/fluid/regularization/bcknd/device/opencl/entropy_viscosity_kernel.cl rename to src/viscous_regularization/avm/bcknd/device/opencl/entropy_viscosity_kernel.cl diff --git a/src/fluid/regularization/entropy_viscosity.f90 b/src/viscous_regularization/avm/entropy_viscosity.f90 similarity index 57% rename from src/fluid/regularization/entropy_viscosity.f90 rename to src/viscous_regularization/avm/entropy_viscosity.f90 index e4868dfd76e8..21c370882a4b 100644 --- a/src/fluid/regularization/entropy_viscosity.f90 +++ b/src/viscous_regularization/avm/entropy_viscosity.f90 @@ -1,4 +1,4 @@ -! Copyright (c) 2025, The Neko Authors +! Copyright (c) 2025-2026, The Neko Authors ! All rights reserved. ! ! Redistribution and use in source and binary forms, with or without @@ -30,9 +30,13 @@ ! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. ! +!> Implements the entropy-based artificial viscosity model. module entropy_viscosity use num_types, only : rp - use regularization, only : regularization_t + use case, only : case_t + use utils, only : neko_error + use registry, only : neko_registry + use artificial_viscosity_model, only : avm_t use json_module, only : json_file use json_utils, only : json_get_or_default use field, only : field_t @@ -40,10 +44,11 @@ module entropy_viscosity use field_series, only : field_series_t use coefs, only : coef_t use dofmap, only : dofmap_t + use gather_scatter, only : gs_t use time_state, only : time_state_t use bdf_time_scheme, only : bdf_time_scheme_t use operators, only : div - use math, only : glmax, absval + use math, only : glmax, absval, col2 use scratch_registry, only : neko_scratch_registry use mesh, only : mesh_t use space, only : space_t @@ -51,7 +56,7 @@ module entropy_viscosity use gs_ops, only : GS_OP_ADD use neko_config, only : NEKO_BCKND_DEVICE use device, only : device_memcpy, HOST_TO_DEVICE, DEVICE_TO_HOST - use device_math, only : device_col3, device_absval, device_glsum + use device_math, only : device_col3, device_absval, device_glsum, device_col2 use entropy_viscosity_cpu, only : entropy_viscosity_compute_residual_cpu, & entropy_viscosity_compute_viscosity_cpu, & entropy_viscosity_apply_element_max_cpu, & @@ -63,83 +68,141 @@ module entropy_viscosity entropy_viscosity_apply_element_max_device, & entropy_viscosity_clamp_to_low_order_device, & entropy_viscosity_smooth_divide_device + use compressible_ops_cpu, only : compressible_ops_cpu_compute_entropy + use compressible_ops_device, only : compressible_ops_device_compute_entropy + use fluid_scheme_compressible, only : fluid_scheme_compressible_t implicit none private - type, public, extends(regularization_t) :: entropy_viscosity_t + !> Entropy-based artificial viscosity model for compressible flow. + type, public, extends(avm_t) :: entropy_viscosity_t + !> Entropy viscosity coefficient. real(kind=rp) :: c_avisc_entropy + !> Low-order viscosity coefficient. real(kind=rp) :: c_avisc_low + !> Ratio of specific heats. + real(kind=rp) :: gamma + !> Entropy residual field. type(field_t) :: entropy_residual + !> Entropy history used by the BDF residual. type(field_series_t) :: S_lag + !> Current entropy field. type(field_t), pointer :: S => null() + !> Pressure field. + type(field_t), pointer :: p => null() + !> Density field. + type(field_t), pointer :: rho => null() + !> First velocity component. type(field_t), pointer :: u => null() + !> Second velocity component. type(field_t), pointer :: v => null() + !> Third velocity component. type(field_t), pointer :: w => null() - type(field_t), pointer :: h => null() + !> Local characteristic mesh size. + type(field_t) :: h + !> Maximum local characteristic wave speed. type(field_t), pointer :: max_wave_speed => null() + !> Computational mesh. type(mesh_t), pointer :: msh => null() + !> Function space. type(space_t), pointer :: Xh => null() + !> Gather-scatter operator for the function space. type(gs_t), pointer :: gs => null() contains + !> Initialize the entropy viscosity model. procedure, pass(this) :: init => entropy_viscosity_init + !> Free the entropy viscosity model. procedure, pass(this) :: free => entropy_viscosity_free - procedure, pass(this) :: compute => entropy_viscosity_compute - procedure, pass(this) :: update_lag => entropy_viscosity_update_lag + !> Compute viscosity before advancing the time step. + procedure, pass(this) :: preprocess => entropy_viscosity_preprocess + !> Advance the entropy history after the time step. + procedure, pass(this) :: compute => entropy_viscosity_update_lag + !> Reconstruct entropy history after restart. + procedure, pass(this) :: restart => entropy_viscosity_restart + !> Compute the local characteristic mesh size. + procedure, pass(this) :: compute_h => entropy_viscosity_compute_h + !> Associate the compressible flow fields. + procedure, pass(this) :: set_fields => entropy_viscosity_set_fields + !> Compute entropy from pressure and density. + procedure, pass(this) :: compute_entropy => & + entropy_viscosity_compute_entropy + !> Compute the entropy residual. procedure, pass(this), private :: compute_residual => & entropy_viscosity_compute_residual + !> Compute the entropy viscosity coefficient. procedure, pass(this), private :: compute_viscosity => & entropy_viscosity_compute_viscosity + !> Smooth viscosity across element interfaces. procedure, pass(this), private :: smooth_viscosity => & entropy_viscosity_smooth_viscosity + !> Replace nodal values by each element maximum. procedure, pass(this), private :: apply_element_max => & entropy_viscosity_apply_element_max + !> Evaluate the low-order viscosity at a point. procedure, pass(this), private :: low_order_viscosity => & entropy_viscosity_low_order end type entropy_viscosity_t - public :: entropy_viscosity_set_fields - contains - subroutine entropy_viscosity_init(this, json, coef, dof, reg_coeff) + !> Initialize the entropy viscosity model from a case and JSON parameters. + !! @param this The entropy viscosity model. + !! @param case The compressible-flow case. + !! @param json The model configuration. + subroutine entropy_viscosity_init(this, case, json) class(entropy_viscosity_t), intent(inout) :: this type(json_file), intent(inout) :: json - type(coef_t), intent(in), target :: coef - type(dofmap_t), intent(in), target :: dof - type(field_t), intent(in), target :: reg_coeff - - call this%init_base(json, coef, dof, reg_coeff) - - call json_get_or_default(json, 'c_avisc_low', this%c_avisc_low, 1.0_rp) + class(case_t), intent(inout), target :: case + character(len=:), allocatable :: reg_coeff_name + + call this%free() + + select type (fluid => case%fluid) + class is (fluid_scheme_compressible_t) + call this%set_fields(fluid%p, fluid%rho, fluid%u, fluid%v, & + fluid%w, fluid%max_wave_speed, fluid%msh, fluid%Xh, fluid%gs_Xh, & + fluid%gamma) + class default + call neko_error('Entropy viscosity requires a compressible fluid scheme') + end select + + call json_get_or_default(json, 'field_name', & + reg_coeff_name, "entropy_viscosity") + call this%init_base(case%fluid%dm_Xh, case%fluid%c_Xh, trim(reg_coeff_name)) + if (allocated(reg_coeff_name)) deallocate(reg_coeff_name) + + call json_get_or_default(json, 'c_avisc_low', this%c_avisc_low, 0.5_rp) call json_get_or_default(json, 'c_avisc_entropy', & this%c_avisc_entropy, 1.0_rp) - call this%entropy_residual%init(dof, 'entropy_residual') + call this%entropy_residual%init(this%dof, 'entropy_residual') - nullify(this%S) - nullify(this%u) - nullify(this%v) - nullify(this%w) - nullify(this%h) - nullify(this%max_wave_speed) - nullify(this%msh) - nullify(this%Xh) - nullify(this%gs) + call neko_registry%add_field(this%dof, 'S', .true.) + this%S => neko_registry%get_field('S') + call this%compute_entropy() + call this%S_lag%init(this%S, 3) + + call this%h%init(this%dof, 'h') + call this%compute_h() end subroutine entropy_viscosity_init + !> Free the entropy viscosity model. + !! @param this The entropy viscosity model. subroutine entropy_viscosity_free(this) class(entropy_viscosity_t), intent(inout) :: this call this%free_base() call this%entropy_residual%free() call this%S_lag%free() + call this%h%free() nullify(this%S) + nullify(this%p) + nullify(this%rho) nullify(this%u) nullify(this%v) nullify(this%w) - nullify(this%h) nullify(this%max_wave_speed) nullify(this%msh) nullify(this%Xh) @@ -147,7 +210,10 @@ subroutine entropy_viscosity_free(this) end subroutine entropy_viscosity_free - subroutine entropy_viscosity_compute(this, time) + !> Compute the artificial viscosity before a time step. + !! @param this The entropy viscosity model. + !! @param time The current time state. + subroutine entropy_viscosity_preprocess(this, time) class(entropy_viscosity_t), intent(inout) :: this type(time_state_t), intent(in) :: time @@ -157,8 +223,13 @@ subroutine entropy_viscosity_compute(this, time) call this%compute_viscosity(time%tstep) - end subroutine entropy_viscosity_compute + end subroutine entropy_viscosity_preprocess + !> Compute the BDF entropy residual. + !! @param this The entropy viscosity model. + !! @param tstep The current time-step index. + !! @param dt The current time-step size. + !! @param dt_lag Previous time-step sizes. subroutine entropy_viscosity_compute_residual(this, tstep, dt, dt_lag) class(entropy_viscosity_t), intent(inout) :: this integer, intent(in) :: tstep @@ -232,6 +303,9 @@ subroutine entropy_viscosity_compute_residual(this, tstep, dt, dt_lag) end subroutine entropy_viscosity_compute_residual + !> Compute and limit the entropy viscosity coefficient. + !! @param this The entropy viscosity model. + !! @param tstep The current time-step index. subroutine entropy_viscosity_compute_viscosity(this, tstep) class(entropy_viscosity_t), intent(inout) :: this integer, intent(in) :: tstep @@ -302,6 +376,7 @@ end subroutine entropy_viscosity_compute_viscosity !> Cross-element smoothing via gather-scatter averaging. !! Averages viscosity values at shared nodes between elements. + !! @param this The entropy viscosity model. subroutine entropy_viscosity_smooth_viscosity(this) class(entropy_viscosity_t), intent(inout) :: this integer :: n @@ -333,6 +408,8 @@ subroutine entropy_viscosity_smooth_viscosity(this) end subroutine entropy_viscosity_smooth_viscosity + !> Replace nodal viscosity values by the maximum in each element. + !! @param this The entropy viscosity model. subroutine entropy_viscosity_apply_element_max(this) class(entropy_viscosity_t), intent(inout) :: this integer :: lx @@ -349,36 +426,92 @@ subroutine entropy_viscosity_apply_element_max(this) end subroutine entropy_viscosity_apply_element_max - subroutine entropy_viscosity_set_fields(this, S, u, v, w, h, max_wave_speed, & - msh, Xh, gs) + !> Associate the fields and discretization used by the model. + !! @param this The entropy viscosity model. + !! @param p The pressure field. + !! @param rho The density field. + !! @param u The first velocity component. + !! @param v The second velocity component. + !! @param w The third velocity component. + !! @param max_wave_speed The maximum characteristic wave speed. + !! @param msh The computational mesh. + !! @param Xh The function space. + !! @param gs The gather-scatter operator. + !! @param gamma The ratio of specific heats. + subroutine entropy_viscosity_set_fields(this, p, rho, u, v, w, & + max_wave_speed, msh, Xh, gs, gamma) class(entropy_viscosity_t), intent(inout) :: this - type(field_t), target, intent(inout) :: S - type(field_t), target, intent(in) :: u, v, w, h, max_wave_speed + type(field_t), target, intent(in) :: p, rho, u, v, w, max_wave_speed type(mesh_t), target, intent(in) :: msh type(space_t), target, intent(in) :: Xh type(gs_t), target, intent(in) :: gs + real(kind=rp), intent(in) :: gamma - this%S => S + this%p => p + this%rho => rho this%u => u this%v => v this%w => w - this%h => h this%max_wave_speed => max_wave_speed this%msh => msh this%Xh => Xh this%gs => gs - - call this%S_lag%init(S, 3) + this%gamma = gamma end subroutine entropy_viscosity_set_fields - subroutine entropy_viscosity_update_lag(this) + !> Compute entropy from the current pressure and density. + !! @param this The entropy viscosity model. + subroutine entropy_viscosity_compute_entropy(this) + class(entropy_viscosity_t), intent(inout) :: this + integer :: n + + n = this%dof%size() + + if (NEKO_BCKND_DEVICE .eq. 1) then + call compressible_ops_device_compute_entropy(this%S, this%p, this%rho, & + this%gamma, n) + else + call compressible_ops_cpu_compute_entropy(this%S%x, this%p%x, & + this%rho%x, this%gamma, n) + end if + + end subroutine entropy_viscosity_compute_entropy + + !> Reinitialize entropy history from a restarted flow state. + !! @param this The entropy viscosity model. + !! @param time The restored time state. + subroutine entropy_viscosity_restart(this, time) + class(entropy_viscosity_t), intent(inout) :: this + type(time_state_t), intent(in) :: time + + call this%compute_entropy() + call this%S_lag%set(this%S) + + end subroutine entropy_viscosity_restart + + !> Shift the entropy history and compute the new entropy. + !! @param this The entropy viscosity model. + !! @param time The current time state. + subroutine entropy_viscosity_update_lag(this, time) class(entropy_viscosity_t), intent(inout) :: this + type(time_state_t), intent(in) :: time call this%S_lag%update() + call this%compute_entropy() + end subroutine entropy_viscosity_update_lag + !> Multiply each velocity component by entropy on the CPU. + !! @param us Entropy-weighted first velocity component. + !! @param vs Entropy-weighted second velocity component. + !! @param ws Entropy-weighted third velocity component. + !! @param u The first velocity component. + !! @param v The second velocity component. + !! @param w The third velocity component. + !! @param S The entropy field. + !! @param n Number of field entries. subroutine entropy_viscosity_col3_vector_cpu(us, vs, ws, u, v, w, S, n) integer, intent(in) :: n real(kind=rp), intent(out) :: us(n), vs(n), ws(n) @@ -398,6 +531,10 @@ subroutine entropy_viscosity_col3_vector_cpu(us, vs, ws, u, v, w, S, n) !$omp end parallel do simd end subroutine entropy_viscosity_col3_vector_cpu + !> Add the flux divergence to the residual and take its absolute value. + !! @param entropy_residual The entropy residual to update. + !! @param div_field The entropy-flux divergence. + !! @param n Number of field entries. subroutine entropy_viscosity_abs_add_cpu(entropy_residual, div_field, n) integer, intent(in) :: n real(kind=rp), intent(inout) :: entropy_residual(n) @@ -415,7 +552,10 @@ subroutine entropy_viscosity_abs_add_cpu(entropy_residual, div_field, n) !$omp end parallel do simd end subroutine entropy_viscosity_abs_add_cpu - !> Compute low-order viscosity at point i: c_avisc_low * h * max_wave_speed + !> Compute low-order viscosity at one point. + !! @param this The entropy viscosity model. + !! @param i The point index. + !! @return The low-order viscosity. pure function entropy_viscosity_low_order(this, i) result(visc) class(entropy_viscosity_t), intent(in) :: this integer, intent(in) :: i @@ -425,4 +565,73 @@ pure function entropy_viscosity_low_order(this, i) result(visc) end function entropy_viscosity_low_order + + !> Compute the characteristic mesh size. + !! Adapted from `les_model_compute_delta` in `les_model.f90`. + !! @todo Move this operation to a shared module. + !! @param this The entropy viscosity model. + subroutine entropy_viscosity_compute_h(this) + class(entropy_viscosity_t), intent(inout) :: this + integer :: e, i, j, k + integer :: im, ip, jm, jp, km, kp + real(kind=rp) :: di, dj, dk, ndim_inv + integer :: lx_half, ly_half, lz_half + + lx_half = this%coef%Xh%lx / 2 + ly_half = this%coef%Xh%ly / 2 + lz_half = this%coef%Xh%lz / 2 + + do concurrent (e = 1:this%coef%msh%nelv) + do concurrent (k = 1:this%coef%Xh%lz, & + j = 1:this%coef%Xh%ly, i = 1:this%coef%Xh%lx) + km = max(1, k-1) + kp = min(this%coef%Xh%lz, k+1) + + jm = max(1, j-1) + jp = min(this%coef%Xh%ly, j+1) + + im = max(1, i-1) + ip = min(this%coef%Xh%lx, i+1) + + di = (this%coef%dof%x(ip, j, k, e) - & + this%coef%dof%x(im, j, k, e))**2 & + + (this%coef%dof%y(ip, j, k, e) - & + this%coef%dof%y(im, j, k, e))**2 & + + (this%coef%dof%z(ip, j, k, e) - & + this%coef%dof%z(im, j, k, e))**2 + + dj = (this%coef%dof%x(i, jp, k, e) - & + this%coef%dof%x(i, jm, k, e))**2 & + + (this%coef%dof%y(i, jp, k, e) - & + this%coef%dof%y(i, jm, k, e))**2 & + + (this%coef%dof%z(i, jp, k, e) - & + this%coef%dof%z(i, jm, k, e))**2 + + dk = (this%coef%dof%x(i, j, kp, e) - & + this%coef%dof%x(i, j, km, e))**2 & + + (this%coef%dof%y(i, j, kp, e) - & + this%coef%dof%y(i, j, km, e))**2 & + + (this%coef%dof%z(i, j, kp, e) - & + this%coef%dof%z(i, j, km, e))**2 + + di = sqrt(di) / (ip - im) + dj = sqrt(dj) / (jp - jm) + dk = sqrt(dk) / (kp - km) + this%h%x(i,j,k,e) = (di * dj * dk)**(1.0_rp / 3.0_rp) + + end do + end do + + if (NEKO_BCKND_DEVICE .eq. 1) then + call device_memcpy(this%h%x, this%h%x_d, this%h%dof%size(),& + HOST_TO_DEVICE, sync = .false.) + call this%gs%op(this%h, GS_OP_ADD) + call device_col2(this%h%x_d, this%coef%mult_d, this%h%dof%size()) + else + call this%gs%op(this%h, GS_OP_ADD) + call col2(this%h%x, this%coef%mult, this%h%dof%size()) + end if + + end subroutine entropy_viscosity_compute_h + end module entropy_viscosity diff --git a/src/viscous_regularization/viscous_regularization.f90 b/src/viscous_regularization/viscous_regularization.f90 new file mode 100644 index 000000000000..6577f3a30fc5 --- /dev/null +++ b/src/viscous_regularization/viscous_regularization.f90 @@ -0,0 +1,153 @@ +! Copyright (c) 2025-2026, The Neko Authors +! All rights reserved. +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! +! * Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! +! * Redistributions in binary form must reproduce the above +! copyright notice, this list of conditions and the following +! disclaimer in the documentation and/or other materials provided +! with the distribution. +! +! * Neither the name of the authors nor the names of its +! contributors may be used to endorse or promote products derived +! from this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE +! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +!> Defines the base interface for viscous regularization methods. +module viscous_regularization + use num_types, only : rp + use json_module, only : json_file + use field, only : field_t + use coefs, only : coef_t + use dofmap, only : dofmap_t + use time_state, only : time_state_t + implicit none + private + + !> Base abstract type for viscous regularization methods. + type, abstract, public :: viscous_regularization_t + !> Spatially varying regularization coefficient. + type(field_t), pointer :: reg_coeff => null() + !> Name of the regularization coefficient field. + character(len=:), allocatable :: reg_coeff_name + !> SEM coefficients. + type(coef_t), pointer :: coef => null() + !> SEM map of degrees of freedom. + type(dofmap_t), pointer :: dof => null() + contains + !> Initialize state shared by all regularization methods. + procedure, pass(this) :: init_base => viscous_regularization_init_base + !> Free state shared by all regularization methods. + procedure, pass(this) :: free_base => viscous_regularization_free_base + !> Initialize a regularization method. + procedure(reg_init), pass(this), deferred :: init + !> Free a regularization method. + procedure(reg_free), pass(this), deferred :: free + !> Update the effective viscosity field. + procedure(reg_update), pass(this), deferred :: update + end type viscous_regularization_t + + abstract interface + !> Initialize a viscous regularization method. + !! @param this The viscous regularization method. + !! @param json The regularization configuration. + !! @param coef The SEM coefficients. + !! @param dof The SEM map of degrees of freedom. + subroutine reg_init(this, json, coef, dof) + import viscous_regularization_t, json_file, coef_t, dofmap_t, field_t + class(viscous_regularization_t), intent(inout) :: this + type(json_file), intent(inout) :: json + type(coef_t), intent(in), target :: coef + type(dofmap_t), intent(in), target :: dof + end subroutine reg_init + end interface + + abstract interface + !> Free a viscous regularization method. + !! @param this The viscous regularization method. + subroutine reg_free(this) + import viscous_regularization_t + class(viscous_regularization_t), intent(inout) :: this + end subroutine reg_free + end interface + + abstract interface + !> Update an effective viscosity field. + !! @param this The viscous regularization method. + !! @param effective_visc The effective viscosity field to update. + !! @param mu The optional physical viscosity field. + subroutine reg_update(this, effective_visc, mu) + import viscous_regularization_t, field_t + class(viscous_regularization_t), intent(inout) :: this + type(field_t), intent(inout) :: effective_visc + type(field_t), intent(in), optional :: mu + end subroutine reg_update + end interface + + interface + !> Allocate and initialize a viscous regularization method. + !! @param object The regularization object to allocate. + !! @param type_name The regularization type name. + !! @param json The regularization configuration. + !! @param coef The SEM coefficients. + !! @param dof The SEM map of degrees of freedom. + module subroutine viscous_regularization_factory(object, type_name, json, & + coef, dof) + class(viscous_regularization_t), allocatable, intent(inout) :: object + character(len=*), intent(in) :: type_name + type(json_file), intent(inout) :: json + type(coef_t), intent(in), target :: coef + type(dofmap_t), intent(in), target :: dof + end subroutine viscous_regularization_factory + end interface + + public :: viscous_regularization_factory + +contains + + !> Initialize state shared by all viscous regularization methods. + !! @param this The viscous regularization method. + !! @param json The regularization configuration. + !! @param coef The SEM coefficients. + !! @param dof The SEM map of degrees of freedom. + subroutine viscous_regularization_init_base(this, json, coef, dof) + class(viscous_regularization_t), intent(inout) :: this + type(json_file), intent(inout) :: json + type(coef_t), intent(in), target :: coef + type(dofmap_t), intent(in), target :: dof + + this%coef => coef + this%dof => dof + + end subroutine viscous_regularization_init_base + + !> Free state shared by all viscous regularization methods. + !! @param this The viscous regularization method. + subroutine viscous_regularization_free_base(this) + class(viscous_regularization_t), intent(inout) :: this + + nullify(this%reg_coeff) + if (allocated(this%reg_coeff_name)) deallocate(this%reg_coeff_name) + nullify(this%coef) + nullify(this%dof) + + end subroutine viscous_regularization_free_base + +end module viscous_regularization diff --git a/src/fluid/regularization_fctry.f90 b/src/viscous_regularization/viscous_regularization_fctry.f90 similarity index 65% rename from src/fluid/regularization_fctry.f90 rename to src/viscous_regularization/viscous_regularization_fctry.f90 index 1481e2d14000..b78efbef1c12 100644 --- a/src/fluid/regularization_fctry.f90 +++ b/src/viscous_regularization/viscous_regularization_fctry.f90 @@ -1,4 +1,4 @@ -! Copyright (c) 2025, The Neko Authors +! Copyright (c) 2025-2026, The Neko Authors ! All rights reserved. ! ! Redistribution and use in source and binary forms, with or without @@ -30,21 +30,27 @@ ! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. ! -submodule(regularization) regularization_fctry - use entropy_viscosity, only : entropy_viscosity_t +!> Implements the viscous regularization factory. +submodule(viscous_regularization) viscous_regularization_fctry + use artificial_viscosity, only : artificial_viscosity_t use utils, only : neko_error implicit none contains - module subroutine regularization_factory(object, type_name, json, & - coef, dof, reg_coeff) - class(regularization_t), allocatable, intent(inout) :: object + !> Allocate and initialize a viscous regularization method. + !! @param object The regularization object to allocate. + !! @param type_name The regularization type name. + !! @param json The regularization configuration. + !! @param coef The SEM coefficients. + !! @param dof The SEM map of degrees of freedom. + module subroutine viscous_regularization_factory(object, type_name, json, & + coef, dof) + class(viscous_regularization_t), allocatable, intent(inout) :: object character(len=*), intent(in) :: type_name type(json_file), intent(inout) :: json type(coef_t), intent(in), target :: coef type(dofmap_t), intent(in), target :: dof - type(field_t), intent(in), target :: reg_coeff if (allocated(object)) then call object%free() @@ -52,14 +58,15 @@ module subroutine regularization_factory(object, type_name, json, & end if select case (trim(type_name)) - case ('entropy', 'entropy_viscosity') - allocate(entropy_viscosity_t::object) + case ('artificial_viscosity') + allocate(artificial_viscosity_t::object) case default - call neko_error('Unknown regularization type: ' // trim(type_name)) + call neko_error('Unknown viscous_regularization type: ' & + // trim(type_name)) end select - call object%init(json, coef, dof, reg_coeff) + call object%init(json, coef, dof) - end subroutine regularization_factory + end subroutine viscous_regularization_factory -end submodule regularization_fctry +end submodule viscous_regularization_fctry diff --git a/tests/integration/tests/test_case_schema/test_case_schema.py b/tests/integration/tests/test_case_schema/test_case_schema.py index c302b0d39ce7..5b3f51aa05b7 100644 --- a/tests/integration/tests/test_case_schema/test_case_schema.py +++ b/tests/integration/tests/test_case_schema/test_case_schema.py @@ -67,6 +67,28 @@ def test_rejects_known_ignored_properties(tmp_path, property_path): assert result.returncode != 0 +@pytest.mark.parametrize( + "property_name", ["c_avisc_low", "c_avisc_entropy"] +) +def test_rejects_legacy_artificial_viscosity_parameters( + tmp_path, property_name +): + with (EXAMPLES_DIR / "tgv" / "tgv.case").open(encoding="utf-8") as handle: + data = json5.load(handle) + + data["case"]["numerics"][property_name] = 0.5 + case_file = tmp_path / "invalid.case" + case_file.write_text(json.dumps(data), encoding="utf-8") + result = subprocess.run( + [sys.executable, str(VALIDATOR), str(case_file)], + capture_output=True, + text=True, + ) + + assert result.returncode != 0 + assert property_name in result.stdout + + def test_allows_user_extension_properties(tmp_path): with (EXAMPLES_DIR / "tgv" / "tgv.case").open(encoding="utf-8") as handle: data = json5.load(handle)