diff --git a/include/LowMachEquationSystem.h b/include/LowMachEquationSystem.h index 8f0282120..29c9ffb7c 100644 --- a/include/LowMachEquationSystem.h +++ b/include/LowMachEquationSystem.h @@ -297,6 +297,7 @@ class ContinuityEquationSystem : public EquationSystem ScalarFieldType* pressure_; VectorFieldType* dpdx_; ScalarFieldType* massFlowRate_; + ScalarFieldType* massVOFBalancedFlowRate_; VectorFieldType* coordinates_; ScalarFieldType* pTmp_; diff --git a/include/SolutionOptions.h b/include/SolutionOptions.h index 6dde276fc..428207316 100644 --- a/include/SolutionOptions.h +++ b/include/SolutionOptions.h @@ -100,6 +100,8 @@ class SolutionOptions bool has_set_boussinesq_time_scale(); + bool realm_has_vof_{false}; + double hybridDefault_; double alphaDefault_; double alphaUpwDefault_; diff --git a/include/VolumeOfFluidEquationSystem.h b/include/VolumeOfFluidEquationSystem.h index ad8e9d2a7..7aaecac98 100644 --- a/include/VolumeOfFluidEquationSystem.h +++ b/include/VolumeOfFluidEquationSystem.h @@ -94,6 +94,8 @@ class VolumeOfFluidEquationSystem : public EquationSystem virtual void solve_and_update(); + virtual void predict_state(); + const bool managePNG_; ScalarFieldType* volumeOfFluid_; VectorFieldType* dvolumeOfFluiddx_; diff --git a/include/edge_kernels/MomentumEdgeSolverAlg.h b/include/edge_kernels/MomentumEdgeSolverAlg.h index c770936ae..61106ef5b 100644 --- a/include/edge_kernels/MomentumEdgeSolverAlg.h +++ b/include/edge_kernels/MomentumEdgeSolverAlg.h @@ -32,7 +32,9 @@ class MomentumEdgeSolverAlg : public AssembleEdgeSolverAlgorithm unsigned dudx_{stk::mesh::InvalidOrdinal}; unsigned edgeAreaVec_{stk::mesh::InvalidOrdinal}; unsigned massFlowRate_{stk::mesh::InvalidOrdinal}; + unsigned massVofBalancedFlowRate_{stk::mesh::InvalidOrdinal}; unsigned viscosity_{stk::mesh::InvalidOrdinal}; + unsigned density_{stk::mesh::InvalidOrdinal}; unsigned pecletFactor_{stk::mesh::InvalidOrdinal}; unsigned maskNodeField_{stk::mesh::InvalidOrdinal}; }; diff --git a/include/edge_kernels/VOFAdvectionEdgeAlg.h b/include/edge_kernels/VOFAdvectionEdgeAlg.h index 79b4dbc4e..d4783de3b 100644 --- a/include/edge_kernels/VOFAdvectionEdgeAlg.h +++ b/include/edge_kernels/VOFAdvectionEdgeAlg.h @@ -39,7 +39,11 @@ class VOFAdvectionEdgeAlg : public AssembleEdgeSolverAlgorithm unsigned dqdx_{stk::mesh::InvalidOrdinal}; unsigned edgeAreaVec_{stk::mesh::InvalidOrdinal}; unsigned massFlowRate_{stk::mesh::InvalidOrdinal}; + unsigned massVofBalancedFlowRate_{stk::mesh::InvalidOrdinal}; unsigned density_{stk::mesh::InvalidOrdinal}; + + double density_liquid_; + double density_gas_; }; } // namespace nalu diff --git a/include/user_functions/DropletVOFAuxFunction.h b/include/user_functions/DropletVOFAuxFunction.h index ba06ccef3..df0e4bc72 100644 --- a/include/user_functions/DropletVOFAuxFunction.h +++ b/include/user_functions/DropletVOFAuxFunction.h @@ -20,7 +20,7 @@ namespace nalu { class DropletVOFAuxFunction : public AuxFunction { public: - DropletVOFAuxFunction(); + DropletVOFAuxFunction(const std::vector& params); virtual ~DropletVOFAuxFunction() {} @@ -34,6 +34,15 @@ class DropletVOFAuxFunction : public AuxFunction const unsigned fieldSize, const unsigned beginPos, const unsigned endPos) const; + + int surf_idx_; + double droppos_x_; + double droppos_y_; + double droppos_z_; + double surf_pos_; + double surf_idx_dbl_; + double radius_; + double interface_thickness_; }; } // namespace nalu diff --git a/include/user_functions/DropletVelocityAuxFunction.h b/include/user_functions/DropletVelocityAuxFunction.h new file mode 100644 index 000000000..ebc26e9dc --- /dev/null +++ b/include/user_functions/DropletVelocityAuxFunction.h @@ -0,0 +1,54 @@ +// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC +// (NTESS), National Renewable Energy Laboratory, University of Texas Austin, +// Northwest Research Associates. Under the terms of Contract DE-NA0003525 +// with NTESS, the U.S. Government retains certain rights in this software. +// +// This software is released under the BSD 3-clause license. See LICENSE file +// for more details. +// + +#ifndef DropletVelocityAuxFunction_h +#define DropletVelocityAuxFunction_h + +#include + +#include + +namespace sierra { +namespace nalu { + +class DropletVelocityAuxFunction : public AuxFunction +{ +public: + DropletVelocityAuxFunction( + const unsigned beginPos, + const unsigned endPos, + const std::vector& params); + + virtual ~DropletVelocityAuxFunction() {} + + using AuxFunction::do_evaluate; + virtual void do_evaluate( + const double* coords, + const double time, + const unsigned spatialDimension, + const unsigned numPoints, + double* fieldPtr, + const unsigned fieldSize, + const unsigned beginPos, + const unsigned endPos) const; + + double droppos_x_; + double droppos_y_; + double droppos_z_; + double dropvel_x_; + double dropvel_y_; + double dropvel_z_; + double radius_; + double interface_thickness_; +}; + +} // namespace nalu +} // namespace sierra + +#endif diff --git a/include/user_functions/SloshingTankVOFAuxFunction.h b/include/user_functions/SloshingTankVOFAuxFunction.h new file mode 100644 index 000000000..f71501de2 --- /dev/null +++ b/include/user_functions/SloshingTankVOFAuxFunction.h @@ -0,0 +1,47 @@ +// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC +// (NTESS), National Renewable Energy Laboratory, University of Texas Austin, +// Northwest Research Associates. Under the terms of Contract DE-NA0003525 +// with NTESS, the U.S. Government retains certain rights in this software. +// +// This software is released under the BSD 3-clause license. See LICENSE file +// for more details. +// + +#ifndef SloshingTankVOFAuxFunction_h +#define SloshingTankVOFAuxFunction_h + +#include + +#include + +namespace sierra { +namespace nalu { + +class SloshingTankVOFAuxFunction : public AuxFunction +{ +public: + SloshingTankVOFAuxFunction(const std::vector& params); + + virtual ~SloshingTankVOFAuxFunction() {} + + using AuxFunction::do_evaluate; + virtual void do_evaluate( + const double* coords, + const double time, + const unsigned spatialDimension, + const unsigned numPoints, + double* fieldPtr, + const unsigned fieldSize, + const unsigned beginPos, + const unsigned endPos) const; + + double water_level_; + double amplitude_; + double kappa_; + double interface_thickness_; +}; + +} // namespace nalu +} // namespace sierra + +#endif diff --git a/reg_tests/CTestList.cmake b/reg_tests/CTestList.cmake index 2ba4247fd..2149800ad 100644 --- a/reg_tests/CTestList.cmake +++ b/reg_tests/CTestList.cmake @@ -193,6 +193,7 @@ if(NOT ENABLE_CUDA AND NOT ENABLE_ROCM) add_test_r(taylorGreenVortex_p3 4) add_test_r(vortexOpen 4) add_test_r(VOFDroplet 4) + add_test_r(VOFInertialDroplet 4) add_test_r(VOFZalDisk 4) if (ENABLE_FFTW) diff --git a/reg_tests/test_files/VOFDroplet/VOFDroplet.norm.gold b/reg_tests/test_files/VOFDroplet/VOFDroplet.norm.gold index 4f2573b88..2b29aadbc 100644 --- a/reg_tests/test_files/VOFDroplet/VOFDroplet.norm.gold +++ b/reg_tests/test_files/VOFDroplet/VOFDroplet.norm.gold @@ -1,3 +1,3 @@ -2.673434136755543e-05 11 0.0001 -6.859619665424398e-06 12 0.0002 -3.313464463487801e-06 13 0.0003 +0.0006181350623480881 11 0.0001 +0.0001148227599440996 12 0.0002 +0.0001728609404847864 13 0.0003 diff --git a/reg_tests/test_files/VOFDroplet/VOFDroplet.yaml b/reg_tests/test_files/VOFDroplet/VOFDroplet.yaml index 6bb2c5226..84905c5ce 100755 --- a/reg_tests/test_files/VOFDroplet/VOFDroplet.yaml +++ b/reg_tests/test_files/VOFDroplet/VOFDroplet.yaml @@ -58,6 +58,8 @@ realms: target_name: block_1 user_function_name: volume_of_fluid: droplet + user_function_parameters: + volume_of_fluid: [0.0, 0.25, 0.0, 0.1, 0.0, 1, 0.0025] - constant: ic_2 target_name: block_1 value: diff --git a/reg_tests/test_files/VOFInertialDroplet/VOFInertialDroplet.norm.gold b/reg_tests/test_files/VOFInertialDroplet/VOFInertialDroplet.norm.gold new file mode 100644 index 000000000..42b655c00 --- /dev/null +++ b/reg_tests/test_files/VOFInertialDroplet/VOFInertialDroplet.norm.gold @@ -0,0 +1,3 @@ +0.0004280828978167022 11 0.005 +0.0002425402745256439 12 0.01 +0.0001207937198159064 13 0.015 diff --git a/reg_tests/test_files/VOFInertialDroplet/VOFInertialDroplet.yaml b/reg_tests/test_files/VOFInertialDroplet/VOFInertialDroplet.yaml new file mode 100755 index 000000000..a2f71f241 --- /dev/null +++ b/reg_tests/test_files/VOFInertialDroplet/VOFInertialDroplet.yaml @@ -0,0 +1,136 @@ +Simulations: + - name: sim1 + time_integrator: ti_1 + optimizer: opt1 + +linear_solvers: + + - name: solve_scalar + type: tpetra + method: gmres + preconditioner: sgs + tolerance: 1e-12 + max_iterations: 50 + kspace: 50 + output_level: 0 + + - name: solve_cont + type: tpetra + method: gmres + preconditioner: muelu + tolerance: 1e-12 + max_iterations: 75 + kspace: 75 + output_level: 0 + recompute_preconditioner: yes + muelu_xml_file_name: ../../xml/vof_resolved.xml + + +realms: + + - name: realm_1 + mesh: "generated:75x75x75|bbox:-0.5,-0.5,-0.5,0.5,0.5,0.5|sideset:xXyYzZ|show" + use_edges: yes + automatic_decomposition_type: rcb + + equation_systems: + name: theEqSys + max_iterations: 3 + + solver_system_specification: + volume_of_fluid: solve_scalar + velocity: solve_scalar + pressure: solve_cont + + systems: + - VolumeOfFluid: + name: myVOF + max_iterations: 1 + convergence_tolerance: 1e-8 + - LowMachEOM: + name: myLowMach + max_iterations: 1 + convergence_tolerance: 1e-8 + + initial_conditions: + + - user_function: ic_1 + target_name: block_1 + user_function_name: + volume_of_fluid: droplet + velocity: droplet + user_function_parameters: + volume_of_fluid: [0.0, 0.0, 0.0, 0.15, -100, 0, 0.01] + velocity: [0.0, 0.0, 0.0, 0.2, 0.2, 0.2, 0.15, 0.01] + - constant: ic_2 + target_name: block_1 + value: + pressure: 0.0 + + material_properties: + target_name: block_1 + specifications: + - name: density + type: volume_of_fluid + primary_value: 1000.0 + secondary_value: 1.0 + + - name: viscosity + type: volume_of_fluid + primary_value: 1.0e-3 + secondary_value: 1.0e-5 + + boundary_conditions: + - periodic_boundary_condition: bc_1 + target_name: [surface_1, surface_2] + periodic_user_data: + search_tolerance: 0.001 + - periodic_boundary_condition: bc_2 + target_name: [surface_3, surface_4] + periodic_user_data: + search_tolerance: 0.001 + - periodic_boundary_condition: bc_3 + target_name: [surface_5, surface_6] + periodic_user_data: + search_tolerance: 0.001 + + solution_options: + name: myOptions + projected_timescale_type: momentum_diag_inv + + options: + - hybrid_factor: + velocity: 0.05 + volume_of_fluid: 0.0 + + - limiter: + volume_of_fluid: yes + pressure: no + velocity: no + + - consistent_mass_matrix_png: + pressure: no + + output: + output_data_base_name: out/inertialDroplet.e + output_frequency: 50 + output_node_set: yes + output_variables: + - density + - volume_of_fluid + - velocity + - pressure + - dvolume_of_fluiddx +Time_Integrators: + - StandardTimeIntegrator: + name: ti_1 + start_time: 0 + termination_time: 0.015 + time_step: 0.005 + time_stepping_type: fixed + time_step_count: 10 + second_order_accuracy: yes + + + realms: + - realm_1 diff --git a/src/LowMachEquationSystem.C b/src/LowMachEquationSystem.C index 28604eb4b..5797d6762 100644 --- a/src/LowMachEquationSystem.C +++ b/src/LowMachEquationSystem.C @@ -174,6 +174,8 @@ #include #include +#include + // deprecated // stk_util @@ -663,6 +665,13 @@ LowMachEquationSystem::register_initial_condition_fcn( theAuxFunc = new SinProfileChannelFlowVelocityAuxFunction(0, nDim); } else if (fcnName == "PerturbedShearLayer") { theAuxFunc = new PerturbedShearLayerVelocityAuxFunction(0, nDim); + } else if (fcnName == "droplet") { + std::map>::const_iterator iterParams = + theParams.find(dofName); + std::vector fcnParams = (iterParams != theParams.end()) + ? (*iterParams).second + : std::vector(); + theAuxFunc = new DropletVelocityAuxFunction(0, nDim, fcnParams); } else { throw std::runtime_error( "InitialCondFunction::non-supported velocity IC"); @@ -2895,6 +2904,12 @@ ContinuityEquationSystem::register_edge_fields( massFlowRate_ = &(meta_data.declare_field( stk::topology::EDGE_RANK, "mass_flow_rate")); stk::mesh::put_field_on_mesh(*massFlowRate_, selector, nullptr); + + if (realm_.solutionOptions_->realm_has_vof_) { + auto massVofBalancedFlowRate_ = &(meta_data.declare_field( + stk::topology::EDGE_RANK, "mass_vof_balanced_flow_rate")); + stk::mesh::put_field_on_mesh(*massVofBalancedFlowRate_, selector, nullptr); + } } //-------------------------------------------------------------------------- diff --git a/src/VolumeOfFluidEquationSystem.C b/src/VolumeOfFluidEquationSystem.C index 846447589..c2ee7c094 100644 --- a/src/VolumeOfFluidEquationSystem.C +++ b/src/VolumeOfFluidEquationSystem.C @@ -7,6 +7,7 @@ // for more details. // +#include #include #include #include @@ -40,14 +41,16 @@ #include "ngp_algorithms/NodalGradEdgeAlg.h" #include "ngp_algorithms/NodalGradBndryElemAlg.h" #include "stk_topology/topology.hpp" -#include "user_functions/ZalesakDiskVOFAuxFunction.h" -#include "user_functions/ZalesakSphereVOFAuxFunction.h" -#include "user_functions/DropletVOFAuxFunction.h" #include "ngp_utils/NgpFieldBLAS.h" #include "ngp_utils/NgpLoopUtils.h" #include "ngp_utils/NgpFieldUtils.h" #include "stk_io/IossBridge.hpp" +#include "user_functions/ZalesakDiskVOFAuxFunction.h" +#include "user_functions/ZalesakSphereVOFAuxFunction.h" +#include "user_functions/DropletVOFAuxFunction.h" +#include "user_functions/SloshingTankVOFAuxFunction.h" + namespace sierra { namespace nalu { @@ -81,6 +84,7 @@ VolumeOfFluidEquationSystem::VolumeOfFluidEquationSystem( // Require div(u) = 0 instead of div(density*u) = 0 realm_.solutionOptions_->solveIncompressibleContinuity_ = true; + realm_.solutionOptions_->realm_has_vof_ = true; // determine nodal gradient form set_nodal_gradient("volume_of_fluid"); @@ -115,7 +119,6 @@ VolumeOfFluidEquationSystem::register_nodal_fields( const int nDim = meta_data.spatial_dimension(); stk::mesh::Selector selector = stk::mesh::selectUnion(part_vec); - // register dof; set it as a restart variable const int numStates = realm_.number_of_states(); auto density_ = &(meta_data.declare_field( @@ -123,7 +126,6 @@ VolumeOfFluidEquationSystem::register_nodal_fields( stk::mesh::put_field_on_mesh(*density_, selector, nullptr); realm_.augment_restart_variable_list("density"); - // push to property list realm_.augment_property_map(DENSITY_ID, density_); volumeOfFluid_ = &(meta_data.declare_field( @@ -131,13 +133,17 @@ VolumeOfFluidEquationSystem::register_nodal_fields( stk::mesh::put_field_on_mesh(*volumeOfFluid_, selector, nullptr); realm_.augment_restart_variable_list("volume_of_fluid"); + auto velocity_ = &(meta_data.declare_field( + stk::topology::NODE_RANK, "velocity", numStates)); + stk::mesh::put_field_on_mesh(*velocity_, selector, nDim, nullptr); + realm_.augment_restart_variable_list("velocity"); + dvolumeOfFluiddx_ = &(meta_data.declare_field( stk::topology::NODE_RANK, "dvolume_of_fluiddx")); stk::mesh::put_field_on_mesh(*dvolumeOfFluiddx_, selector, nDim, nullptr); stk::io::set_field_output_type( *dvolumeOfFluiddx_, stk::io::FieldOutputType::VECTOR_3D); - // delta solution for linear solver; share delta since this is a split system vofTmp_ = &(meta_data.declare_field(stk::topology::NODE_RANK, "vofTmp")); stk::mesh::put_field_on_mesh(*vofTmp_, selector, nullptr); @@ -178,6 +184,9 @@ VolumeOfFluidEquationSystem::register_edge_fields( auto massFlowRate_ = &(meta_data.declare_field( stk::topology::EDGE_RANK, "mass_flow_rate")); stk::mesh::put_field_on_mesh(*massFlowRate_, selector, nullptr); + auto massVofBalancedFlowRate_ = &(meta_data.declare_field( + stk::topology::EDGE_RANK, "mass_vof_balanced_flow_rate")); + stk::mesh::put_field_on_mesh(*massVofBalancedFlowRate_, selector, nullptr); } //-------------------------------------------------------------------------- @@ -262,7 +271,6 @@ VolumeOfFluidEquationSystem::register_inflow_bc( const stk::topology& /* partTopo */, const InflowBoundaryConditionData& inflowBCData) { - // algorithm type const AlgorithmType algType = INFLOW; ScalarFieldType& vofNp1 = volumeOfFluid_->field_of_state(stk::mesh::StateNP1); VectorFieldType& dvofdxNone = @@ -270,12 +278,10 @@ VolumeOfFluidEquationSystem::register_inflow_bc( stk::mesh::MetaData& meta_data = realm_.meta_data(); - // register boundary data; gamma_bc ScalarFieldType* theBcField = &(meta_data.declare_field(stk::topology::NODE_RANK, "vof_bc")); stk::mesh::put_field_on_mesh(*theBcField, *part, nullptr); - // extract the value for user specified tke and save off the AuxFunction InflowUserData userData = inflowBCData.userData_; std::string vofName = "volume_of_fluid"; UserDataType theDataType = get_bc_data_type(userData, vofName); @@ -296,29 +302,22 @@ VolumeOfFluidEquationSystem::register_inflow_bc( "only constant functions supported"); } - // bc data alg AuxFunctionAlgorithm* auxAlg = new AuxFunctionAlgorithm( realm_, part, theBcField, theAuxFunc, stk::topology::NODE_RANK); - // how to populate the field? if (userData.externalData_) { - // xfer will handle population; only need to populate the initial value realm_.initCondAlg_.push_back(auxAlg); } else { - // put it on bcData bcDataAlg_.push_back(auxAlg); } - // copy vof_bc to gamma_transition np1... CopyFieldAlgorithm* theCopyAlg = new CopyFieldAlgorithm( realm_, part, theBcField, &vofNp1, 0, 1, stk::topology::NODE_RANK); bcDataMapAlg_.push_back(theCopyAlg); - // non-solver; dgamdx; allow for element-based shifted nodalGradAlgDriver_.register_face_algorithm( algType, part, "vof_nodal_grad", &vofNp1, &dvofdxNone, edgeNodalGradient_); - // Dirichlet bc std::map::iterator itd = solverAlgDriver_->solverDirichAlgMap_.find(algType); if (itd == solverAlgDriver_->solverDirichAlgMap_.end()) { @@ -347,8 +346,7 @@ VolumeOfFluidEquationSystem::register_open_bc( // non-solver; dvofdx; allow for element-based shifted nodalGradAlgDriver_.register_face_algorithm( - algType, part, "gamma_nodal_grad", &vofNp1, &dvofdxNone, - edgeNodalGradient_); + algType, part, "vof_nodal_grad", &vofNp1, &dvofdxNone, edgeNodalGradient_); } //-------------------------------------------------------------------------- @@ -358,10 +356,9 @@ void VolumeOfFluidEquationSystem::register_wall_bc( stk::mesh::Part* part, const stk::topology& /*theTopo*/, - const WallBoundaryConditionData& /* wallBCData */) + const WallBoundaryConditionData& /*wallBCData*/) { - // algorithm type - const AlgorithmType algType = WALL; + const AlgorithmType algType = SYMMETRY; ScalarFieldType& vofNp1 = volumeOfFluid_->field_of_state(stk::mesh::StateNP1); VectorFieldType& dvofdxNone = @@ -369,8 +366,7 @@ VolumeOfFluidEquationSystem::register_wall_bc( // non-solver; dvofdx; allow for element-based shifted nodalGradAlgDriver_.register_face_algorithm( - algType, part, "gamma_nodal_grad", &vofNp1, &dvofdxNone, - edgeNodalGradient_); + algType, part, "vof_nodal_grad", &vofNp1, &dvofdxNone, edgeNodalGradient_); } //-------------------------------------------------------------------------- @@ -382,7 +378,6 @@ VolumeOfFluidEquationSystem::register_symmetry_bc( const stk::topology& /*theTopo*/, const SymmetryBoundaryConditionData& /* symmetryBCData */) { - // algorithm type const AlgorithmType algType = SYMMETRY; ScalarFieldType& vofNp1 = volumeOfFluid_->field_of_state(stk::mesh::StateNP1); @@ -391,8 +386,7 @@ VolumeOfFluidEquationSystem::register_symmetry_bc( // non-solver; dvofdx; allow for element-based shifted nodalGradAlgDriver_.register_face_algorithm( - algType, part, "gamma_nodal_grad", &vofNp1, &dvofdxNone, - edgeNodalGradient_); + algType, part, "vof_nodal_grad", &vofNp1, &dvofdxNone, edgeNodalGradient_); } //-------------------------------------------------------------------------- @@ -469,7 +463,7 @@ void VolumeOfFluidEquationSystem::register_initial_condition_fcn( stk::mesh::Part* part, const std::map& theNames, - const std::map>& /* theParams */) + const std::map>& theParams) { // iterate map and check for name const std::string dofName = "volume_of_fluid"; @@ -487,14 +481,20 @@ VolumeOfFluidEquationSystem::register_initial_condition_fcn( TurbulenceModel::SST_AMS) ? true : false; + VectorFieldType* velocity_ = realm_.meta_data().get_field( + stk::topology::NODE_RANK, "velocity"); ScalarFieldType* density_ = realm_.meta_data().get_field( stk::topology::NODE_RANK, "density"); - std::vector userSpec(1); - userSpec[0] = 1.0; + + std::vector userSpec = {1.0, 0.0, 0.0}; AuxFunction* constantAuxFunc = new ConstantAuxFunction(0, 1, userSpec); AuxFunctionAlgorithm* constantAuxAlg = new AuxFunctionAlgorithm( realm_, part, density_, constantAuxFunc, stk::topology::NODE_RANK); realm_.initCondAlg_.push_back(constantAuxAlg); + AuxFunction* constantVelFunc = new ConstantAuxFunction(0, 3, userSpec); + AuxFunctionAlgorithm* constantVelAlg = new AuxFunctionAlgorithm( + realm_, part, velocity_, constantVelFunc, stk::topology::NODE_RANK); + realm_.initCondAlg_.push_back(constantVelAlg); auto VOFSetMassFlowRate = new ZalesakDiskMassFlowRateEdgeAlg(realm_, part, this, useAvgMdot); realm_.initCondAlg_.push_back(VOFSetMassFlowRate); @@ -520,7 +520,19 @@ VolumeOfFluidEquationSystem::register_initial_condition_fcn( realm_.initCondAlg_.push_back(VOFSetMassFlowRate); } } else if (fcnName == "droplet") { - theAuxFunc = new DropletVOFAuxFunction(); + std::map>::const_iterator iterParams = + theParams.find(dofName); + std::vector fcnParams = (iterParams != theParams.end()) + ? (*iterParams).second + : std::vector(); + theAuxFunc = new DropletVOFAuxFunction(fcnParams); + } else if (fcnName == "sloshing_tank") { + std::map>::const_iterator iterParams = + theParams.find(dofName); + std::vector fcnParams = (iterParams != theParams.end()) + ? (*iterParams).second + : std::vector(); + theAuxFunc = new SloshingTankVOFAuxFunction(fcnParams); } else { throw std::runtime_error("VolumeOfFluidEquationSystem::register_initial_" "condition_fcn: limited functions supported"); @@ -588,6 +600,7 @@ VolumeOfFluidEquationSystem::compute_projected_nodal_gradient() projectedNodalGradEqs_->solve_and_update_external(); } } + //-------------------------------------------------------------------------- //-------- solve_and_update ------------------------------------------------ //-------------------------------------------------------------------------- @@ -606,7 +619,6 @@ VolumeOfFluidEquationSystem::solve_and_update() NaluEnv::self().naluOutputP0() << " " << k + 1 << "/" << maxIterations_ << std::setw(15) << std::right << userSuppliedName_ << std::endl; - assemble_and_solve(vofTmp_); solution_update(1.0, *vofTmp_, 1.0, *volumeOfFluid_); @@ -614,5 +626,32 @@ VolumeOfFluidEquationSystem::solve_and_update() } } +//-------------------------------------------------------------------------- +//-------- predict_state --------------------------------------------------- +//-------------------------------------------------------------------------- +void +VolumeOfFluidEquationSystem::predict_state() +{ + // const auto& meshInfo = realm_.mesh_info(); + const auto& ngpMesh = realm_.ngp_mesh(); + const auto& fieldMgr = realm_.ngp_field_manager(); + + auto& vofN = fieldMgr.get_field( + volumeOfFluid_->field_of_state(stk::mesh::StateN).mesh_meta_data_ordinal()); + auto& vofNp1 = fieldMgr.get_field( + volumeOfFluid_->field_of_state(stk::mesh::StateNP1) + .mesh_meta_data_ordinal()); + + vofN.sync_to_device(); + + const auto& meta = realm_.meta_data(); + const stk::mesh::Selector sel = + (meta.locally_owned_part() | meta.globally_shared_part() | + meta.aura_part()) & + stk::mesh::selectField(*volumeOfFluid_); + nalu_ngp::field_copy(ngpMesh, sel, vofNp1, vofN, 1); + vofNp1.modify_on_device(); +} + } // namespace nalu } // namespace sierra diff --git a/src/edge_kernels/ContinuityEdgeSolverAlg.C b/src/edge_kernels/ContinuityEdgeSolverAlg.C index 5ce8e31d1..11b777835 100644 --- a/src/edge_kernels/ContinuityEdgeSolverAlg.C +++ b/src/edge_kernels/ContinuityEdgeSolverAlg.C @@ -11,6 +11,7 @@ #include "utils/StkHelpers.h" #include "stk_mesh/base/NgpField.hpp" #include "stk_mesh/base/Types.hpp" +#include "SolutionOptions.h" namespace sierra { namespace nalu { @@ -89,6 +90,7 @@ ContinuityEdgeSolverAlg::execute() const stk::mesh::FastMeshIndex& nodeR) { // Scratch work array for edgeAreaVector NALU_ALIGNED DblType av[NDimMax_]; + // Populate area vector work array for (int d = 0; d < ndim; ++d) av[d] = edgeAreaVec.get(edge, d); @@ -131,7 +133,7 @@ ContinuityEdgeSolverAlg::execute() const DblType ujIp = 0.5 * (velocity.get(nodeR, d) + velocity.get(nodeL, d)); const DblType GjIp = - 0.5 * (Gpdx.get(nodeR, d) / udiagR + Gpdx.get(nodeL, d) / udiagL); + 0.5 * (Gpdx.get(nodeR, d) / (udiagR) + Gpdx.get(nodeL, d) / (udiagL)); tmdot += (interpTogether * rhoUjIp + om_interpTogether * rhoIp * ujIp + GjIp) * av[d] - diff --git a/src/edge_kernels/MomentumEdgeSolverAlg.C b/src/edge_kernels/MomentumEdgeSolverAlg.C index d58293603..2ea4532a5 100644 --- a/src/edge_kernels/MomentumEdgeSolverAlg.C +++ b/src/edge_kernels/MomentumEdgeSolverAlg.C @@ -7,6 +7,7 @@ // for more details. // +#include #include "edge_kernels/MomentumEdgeSolverAlg.h" #include "EquationSystem.h" #include "PecletFunction.h" @@ -41,16 +42,29 @@ MomentumEdgeSolverAlg::MomentumEdgeSolverAlg( else viscName = "viscosity"; + density_ = get_field_ordinal(meta, "density", stk::mesh::StateNP1); viscosity_ = get_field_ordinal(meta, viscName); dudx_ = get_field_ordinal(meta, "dudx"); edgeAreaVec_ = get_field_ordinal(meta, "edge_area_vector", stk::topology::EDGE_RANK); massFlowRate_ = get_field_ordinal(meta, "mass_flow_rate", stk::topology::EDGE_RANK); + massVofBalancedFlowRate_ = + realm.solutionOptions_->realm_has_vof_ + ? get_field_ordinal( + meta, "mass_vof_balanced_flow_rate", stk::topology::EDGE_RANK) + : massFlowRate_; pecletFactor_ = get_field_ordinal(meta, "peclet_factor", stk::topology::EDGE_RANK); maskNodeField_ = get_field_ordinal( meta, "abl_wall_no_slip_wall_func_node_mask", stk::topology::NODE_RANK); + + if (realm_.solutionOptions_->realm_has_vof_) { + NaluEnv::self().naluOutputP0() + << "WARNING: volume_of_fluid is present. For stability, upwinding in the " + "MomentumEdgeSolverAlg is automatically turned on near the liquid-gas " + "interface.\n"; + } } void @@ -62,14 +76,16 @@ MomentumEdgeSolverAlg::execute() const std::string dofName = "velocity"; const DblType includeDivU = realm_.get_divU(); const DblType alpha = realm_.get_alpha_factor(dofName); - const DblType alphaUpw = realm_.get_alpha_upw_factor(dofName); + const DblType alphaUpw_input = realm_.get_alpha_upw_factor(dofName); const DblType hoUpwind = realm_.get_upw_factor(dofName); const DblType relaxFacU = realm_.solutionOptions_->get_relaxation_factor(dofName); const bool useLimiter = realm_.primitive_uses_limiter(dofName); const DblType om_alpha = 1.0 - alpha; - const DblType om_alphaUpw = 1.0 - alphaUpw; + const DblType om_alphaUpw_input = 1.0 - alphaUpw_input; + + const DblType has_vof = (double)realm_.solutionOptions_->realm_has_vof_; // STK stk::mesh::NgpField instances for capture by lambda const auto& fieldMgr = realm_.ngp_field_manager(); @@ -78,8 +94,11 @@ MomentumEdgeSolverAlg::execute() const auto vel = fieldMgr.get_field(velocity_); const auto dudx = fieldMgr.get_field(dudx_); const auto viscosity = fieldMgr.get_field(viscosity_); + const auto density = fieldMgr.get_field(density_); const auto edgeAreaVec = fieldMgr.get_field(edgeAreaVec_); const auto massFlowRate = fieldMgr.get_field(massFlowRate_); + const auto massVofBalancedFlowRate = + fieldMgr.get_field(massVofBalancedFlowRate_); const auto pecletFactor = fieldMgr.get_field(pecletFactor_); const auto maskNodeField = fieldMgr.get_field(maskNodeField_); @@ -89,14 +108,23 @@ MomentumEdgeSolverAlg::execute() ShmemDataType & smdata, const stk::mesh::FastMeshIndex& edge, const stk::mesh::FastMeshIndex& nodeL, const stk::mesh::FastMeshIndex& nodeR) { + // Variables are read-only when C++ captures by value into a lambda + // This avoids this issue, allowing modification of upwinding parameters + // for VOF + DblType alphaUpw = alphaUpw_input; + DblType om_alphaUpw = om_alphaUpw_input; + // Scratch work array for edgeAreaVector NALU_ALIGNED DblType av[NDimMax_]; // Populate area vector work array for (int d = 0; d < ndim; ++d) av[d] = edgeAreaVec.get(edge, d); - const DblType mdot = massFlowRate.get(edge, 0); + const DblType mdot = + massFlowRate.get(edge, 0) + has_vof * massVofBalancedFlowRate(edge, 0); + const DblType densityL = density.get(nodeL, 0); + const DblType densityR = density.get(nodeR, 0); const DblType viscosityL = viscosity.get(nodeL, 0); const DblType viscosityR = viscosity.get(nodeR, 0); @@ -130,9 +158,6 @@ MomentumEdgeSolverAlg::execute() } } - const DblType pecfac = pecletFactor.get(edge, 0); - const DblType om_pecfac = 1.0 - pecfac; - NALU_ALIGNED DblType limitL[NDimMax_] = {1.0, 1.0, 1.0}; NALU_ALIGNED DblType limitR[NDimMax_] = {1.0, 1.0, 1.0}; @@ -146,12 +171,35 @@ MomentumEdgeSolverAlg::execute() } } + // Upwinding switch for multiphase cases: + // Factors determined by ensuring full upwinding at interfaces with + // interface widths that are ~2 cells thick + DblType pecfac = pecletFactor.get(edge, 0); + DblType om_pecfac = 1.0 - pecfac; + DblType density_upwinding_factor = 1.0; + if (has_vof > 0.5) { + const DblType min_density = stk::math::min(densityL, densityR); + const DblType density_differential = + stk::math::abs(densityL - densityR) / min_density; + density_upwinding_factor = + 1.0 - stk::math::erf(6.0 * density_differential); + + alphaUpw = density_upwinding_factor * alphaUpw + + (1.0 - density_upwinding_factor); + om_alphaUpw = 1.0 - alphaUpw; + pecfac = + 1.0 - density_upwinding_factor + density_upwinding_factor * pecfac; + om_pecfac = 1.0 - pecfac; + } + // Upwind extrapolation with limiter terms NALU_ALIGNED DblType uIpL[NDimMax_]; NALU_ALIGNED DblType uIpR[NDimMax_]; for (int d = 0; d < ndim; ++d) { - uIpL[d] = vel.get(nodeL, d) + duL[d] * hoUpwind * limitL[d]; - uIpR[d] = vel.get(nodeR, d) - duR[d] * hoUpwind * limitR[d]; + uIpL[d] = vel.get(nodeL, d) + + duL[d] * hoUpwind * limitL[d] * density_upwinding_factor; + uIpR[d] = vel.get(nodeR, d) - + duR[d] * hoUpwind * limitR[d] * density_upwinding_factor; } // TODO(psakiev) extract this a funciton into EdgeKernelUtils.h diff --git a/src/edge_kernels/VOFAdvectionEdgeAlg.C b/src/edge_kernels/VOFAdvectionEdgeAlg.C index 2a49c85fb..0b07f112d 100644 --- a/src/edge_kernels/VOFAdvectionEdgeAlg.C +++ b/src/edge_kernels/VOFAdvectionEdgeAlg.C @@ -16,6 +16,9 @@ #include "stk_mesh/base/NgpField.hpp" #include "stk_mesh/base/Types.hpp" #include +#include +#include +#include "ngp_utils/NgpLoopUtils.h" namespace sierra { namespace nalu { @@ -39,19 +42,50 @@ VOFAdvectionEdgeAlg::VOFAdvectionEdgeAlg( massFlowRate_ = get_field_ordinal( meta, (useAverages) ? "average_mass_flow_rate" : "mass_flow_rate", stk::topology::EDGE_RANK); + massVofBalancedFlowRate_ = get_field_ordinal( + meta, "mass_vof_balanced_flow_rate", stk::topology::EDGE_RANK); density_ = get_field_ordinal(realm.meta_data(), "density", stk::mesh::StateNP1); + + std::map::iterator itf = + realm_.materialPropertys_.propertyDataMap_.find(DENSITY_ID); + + // Hard set value here for unit testing without property map. + if (itf == realm_.materialPropertys_.propertyDataMap_.end()) { + density_liquid_ = 1000.0; + density_gas_ = 1.0; + } else { + auto mdata = (*itf).second; + + switch (mdata->type_) { + case CONSTANT_MAT: { + density_liquid_ = mdata->constValue_; + density_gas_ = mdata->constValue_; + break; + } + case VOF_MAT: { + density_liquid_ = mdata->primary_; + density_gas_ = mdata->secondary_; + break; + } + default: { + throw std::runtime_error("Incorrect density property set for VOF " + "calculations. Use a constant or " + "VOF property for density."); + break; + } + } + } } void VOFAdvectionEdgeAlg::execute() { - const double eps = 1.0e-16; + const double eps = 1.0e-11; const double gradient_eps = 1.0e-9; - // Could be made into user paramter for more control. - const double compression_magnitude = 0.1; const int ndim = realm_.meta_data().spatial_dimension(); + const auto& meta = realm_.meta_data(); const DblType alphaUpw = realm_.get_alpha_upw_factor("volume_of_fluid"); const DblType hoUpwind = realm_.get_upw_factor("volume_of_fluid"); @@ -68,6 +102,8 @@ VOFAdvectionEdgeAlg::execute() const auto dqdx = fieldMgr.get_field(dqdx_); const auto edgeAreaVec = fieldMgr.get_field(edgeAreaVec_); const auto massFlowRate = fieldMgr.get_field(massFlowRate_); + const auto massVofBalancedFlowRate = + fieldMgr.get_field(massVofBalancedFlowRate_); const auto density = fieldMgr.get_field(density_); run_algorithm( @@ -125,7 +161,6 @@ VOFAdvectionEdgeAlg::execute() // Advective flux const DblType qIp = 0.5 * (qNp1R + qNp1L); // 2nd order central term - // Upwinded term const DblType qUpw = (vdot > 0) ? (alphaUpw * qIpL + om_alphaUpw * qIp) : (alphaUpw * qIpR + om_alphaUpw * qIp); @@ -145,47 +180,118 @@ VOFAdvectionEdgeAlg::execute() smdata.lhs(0, 1) += alhsfac; // Compression term - const DblType velocity_scale = stk::math::abs( - vdot / stk::math::sqrt(av[0] * av[0] + av[1] * av[1] + av[2] * av[2])); - DblType dqdxMagL = 0.0; - DblType dqdxMagR = 0.0; + DblType dOmegadxMag = 0.0; + DblType interface_gradient[3] = {0.0, 0.0, 0.0}; + + for (int j = 0; j < ndim; ++j) { + interface_gradient[j] = 0.5 * (dqdx.get(nodeL, j) + dqdx.get(nodeR, j)); + } DblType interface_normal[3] = {0.0, 0.0, 0.0}; + for (int j = 0; j < ndim; ++j) + dOmegadxMag += interface_gradient[j] * interface_gradient[j]; + + dOmegadxMag = stk::math::sqrt(dOmegadxMag); + + // No gradient == no interface + if (dOmegadxMag < gradient_eps) + return; + + for (int d = 0; d < ndim; ++d) + interface_normal[d] = interface_gradient[d] / dOmegadxMag; + + DblType axdx = 0.0; + DblType asq = 0.0; + DblType diffusion_coef = 0.0; + + for (int d = 0; d < ndim; ++d) { + const DblType dxj = + coordinates.get(nodeR, d) - coordinates.get(nodeL, d); + diffusion_coef += dxj * dxj; + asq += av[d] * av[d]; + axdx += av[d] * dxj; + } + + // Hard-coded values comes from Jain, 2022 to enforce + // VOF function bounds of [0,1] while maintaining interface + // thickness that is ~2 cells + + const DblType velocity_scale = + 5.0 * stk::math::abs( + vdot / + stk::math::sqrt(av[0] * av[0] + av[1] * av[1] + av[2] * av[2])); + + diffusion_coef = stk::math::sqrt(diffusion_coef) * 0.3; + + const DblType inv_axdx = 1.0 / axdx; + + const DblType dlhsfac = -velocity_scale * diffusion_coef * asq * inv_axdx; + + smdata.rhs(0) -= dlhsfac * (qNp1R - qNp1L); + smdata.rhs(1) += dlhsfac * (qNp1R - qNp1L); + + massVofBalancedFlowRate.get(edge, 0) = + dlhsfac * (qNp1R - qNp1L) * (density_liquid_ - density_gas_); + + smdata.lhs(0, 0) -= dlhsfac; + smdata.lhs(0, 1) += dlhsfac; + + smdata.lhs(1, 0) += dlhsfac; + smdata.lhs(1, 1) -= dlhsfac; + + const DblType omegaL = + diffusion_coef * stk::math::log((qNp1L + eps) / (1.0 - qNp1L + eps)); + const DblType omegaR = + diffusion_coef * stk::math::log((qNp1R + eps) / (1.0 - qNp1R + eps)); + const DblType omegaIp = 0.5 * (omegaL + omegaR); + + dOmegadxMag = 0.0; + + for (int d = 0; d < 3; ++d) { + interface_gradient[d] = 0.0; + interface_normal[d] = 0.0; + } + for (int j = 0; j < ndim; ++j) { - dqdxMagL += dqdx.get(nodeL, j) * dqdx.get(nodeL, j); - dqdxMagR += dqdx.get(nodeR, j) * dqdx.get(nodeR, j); - interface_normal[j] = - 0.5 * dqdx.get(nodeL, j) + 0.5 * dqdx.get(nodeR, j); + interface_gradient[j] = 0.5 * (dqdx.get(nodeL, j) + dqdx.get(nodeR, j)); + interface_gradient[j] *= (2.0 * diffusion_coef * eps + diffusion_coef) / + (eps * eps + eps - qIp * qIp + qIp); } - dqdxMagL = stk::math::sqrt(dqdxMagL); - dqdxMagR = stk::math::sqrt(dqdxMagR); + for (int j = 0; j < ndim; ++j) + dOmegadxMag += interface_gradient[j] * interface_gradient[j]; + + dOmegadxMag = stk::math::sqrt(dOmegadxMag); // No gradient == no interface - if ( - stk::math::abs(dqdxMagL) + stk::math::abs(dqdxMagR) < - 2.0 * gradient_eps) + if (dOmegadxMag < gradient_eps) return; for (int d = 0; d < ndim; ++d) - interface_normal[d] /= 0.5 * dqdxMagL + 0.5 * dqdxMagR + eps; + interface_normal[d] = interface_gradient[d] / dOmegadxMag; DblType compression = 0.0; for (int d = 0; d < ndim; ++d) - compression += compression_magnitude * interface_normal[d] * - velocity_scale * qIp * (1.0 - qIp) * av[d]; + compression += + velocity_scale * 0.25 * + (1.0 - stk::math::tanh(0.5 * omegaIp / diffusion_coef) * + stk::math::tanh(0.5 * omegaIp / diffusion_coef)) * + interface_normal[d] * av[d]; smdata.rhs(0) -= compression; smdata.rhs(1) += compression; + massVofBalancedFlowRate.get(edge, 0) += + compression * (density_liquid_ - density_gas_); + // Left node contribution; Lag in iterations except for central 0.5*q term DblType slhsfac = 0.0; for (int d = 0; d < ndim; ++d) - slhsfac += compression_magnitude * 0.5 * interface_normal[d] * - velocity_scale * (1.0 - qIp) * av[d]; + slhsfac += 0.5 * interface_normal[d] * 1.5 * velocity_scale * + (1.0 - qIp) * av[d]; smdata.lhs(0, 0) += slhsfac / relaxFac; smdata.lhs(1, 0) -= slhsfac; diff --git a/src/ngp_algorithms/MdotEdgeAlg.C b/src/ngp_algorithms/MdotEdgeAlg.C index b8e60b7b6..02a454a2e 100644 --- a/src/ngp_algorithms/MdotEdgeAlg.C +++ b/src/ngp_algorithms/MdotEdgeAlg.C @@ -12,6 +12,7 @@ #include "ngp_utils/NgpFieldOps.h" #include "ngp_utils/NgpFieldManager.h" #include "Realm.h" +#include "SolutionOptions.h" #include "utils/StkHelpers.h" #include "stk_mesh/base/NgpMesh.hpp" #include "stk_mesh/base/NgpField.hpp" diff --git a/src/user_functions/CMakeLists.txt b/src/user_functions/CMakeLists.txt index 331624d18..2be2556c2 100644 --- a/src/user_functions/CMakeLists.txt +++ b/src/user_functions/CMakeLists.txt @@ -42,4 +42,6 @@ target_sources(nalu PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/ZalesakSphereVOFAuxFunction.C ${CMAKE_CURRENT_SOURCE_DIR}/ZalesakSphereMassFlowRateKernel.C ${CMAKE_CURRENT_SOURCE_DIR}/DropletVOFAuxFunction.C + ${CMAKE_CURRENT_SOURCE_DIR}/DropletVelocityAuxFunction.C + ${CMAKE_CURRENT_SOURCE_DIR}/SloshingTankVOFAuxFunction.C ) diff --git a/src/user_functions/DropletVOFAuxFunction.C b/src/user_functions/DropletVOFAuxFunction.C index dd075daec..90a8e5c19 100644 --- a/src/user_functions/DropletVOFAuxFunction.C +++ b/src/user_functions/DropletVOFAuxFunction.C @@ -18,9 +18,39 @@ namespace sierra { namespace nalu { -DropletVOFAuxFunction::DropletVOFAuxFunction() : AuxFunction(0, 1) +DropletVOFAuxFunction::DropletVOFAuxFunction(const std::vector& params) + : AuxFunction(0, 1), + surf_idx_(1), + droppos_x_(0.0), + droppos_y_(0.0), + droppos_z_(0.0), + radius_(0.1), + surf_pos_(0.0), + surf_idx_dbl_(1.0), + interface_thickness_(0.0025) { - // does nothing + // check size and populate + if (params.size() != 7 && !params.empty()) + throw std::runtime_error("Realm::setup_initial_conditions: " + "droplet (volume_of_fluid) requires 7 params: 3 " + "components of droplet position, droplet " + "radius, surface position, surface coordinate " + "index, and interface thickness"); + if (!params.empty()) { + droppos_x_ = params[0]; + droppos_y_ = params[1]; + droppos_z_ = params[2]; + radius_ = params[3]; + surf_pos_ = params[4]; + surf_idx_dbl_ = params[5]; + interface_thickness_ = params[6]; + } + + if (surf_idx_dbl_ < 0.5) { + surf_idx_ = 0; + } else if (surf_idx_dbl_ > 1.5) { + surf_idx_ = 2; + } } void @@ -39,13 +69,23 @@ DropletVOFAuxFunction::do_evaluate( const double x = coords[0]; const double y = coords[1]; const double z = coords[2]; - const double interface_thickness = 0.025; - fieldPtr[0] = 0.0; - fieldPtr[0] += -0.5 * (std::erf(y / interface_thickness) + 1.0) + 1.0; + // Flat surface + fieldPtr[0] = + -0.5 * (std::erf((coords[surf_idx_] - surf_pos_) / interface_thickness_) + + 1.0) + + 1.0; + + // Droplet + auto rad_pos = std::sqrt( + (x - droppos_x_) * (x - droppos_x_) + + (y - droppos_y_) * (y - droppos_y_) + + (z - droppos_z_) * (z - droppos_z_)) - + radius_; + fieldPtr[0] += + -0.5 * (std::erf(rad_pos / interface_thickness_) + 1.0) + 1.0; - auto radius = std::sqrt(x * x + (y - 0.25) * (y - 0.25) + z * z) - 0.075; - fieldPtr[0] += -0.5 * (std::erf(radius / interface_thickness) + 1.0) + 1.0; + fieldPtr[0] = std::max(0.0, std::min(1.0, fieldPtr[0])); fieldPtr += fieldSize; coords += spatialDimension; diff --git a/src/user_functions/DropletVelocityAuxFunction.C b/src/user_functions/DropletVelocityAuxFunction.C new file mode 100644 index 000000000..2f92b75ce --- /dev/null +++ b/src/user_functions/DropletVelocityAuxFunction.C @@ -0,0 +1,89 @@ +// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC +// (NTESS), National Renewable Energy Laboratory, University of Texas Austin, +// Northwest Research Associates. Under the terms of Contract DE-NA0003525 +// with NTESS, the U.S. Government retains certain rights in this software. +// +// This software is released under the BSD 3-clause license. See LICENSE file +// for more details. +// + +#include +#include + +// basic c++ +#include +#include +#include + +namespace sierra { +namespace nalu { + +DropletVelocityAuxFunction::DropletVelocityAuxFunction( + const unsigned beginPos, + const unsigned endPos, + const std::vector& params) + : AuxFunction(beginPos, endPos), + droppos_x_(0.0), + droppos_y_(0.0), + droppos_z_(0.0), + dropvel_x_(0.1), + dropvel_y_(0.1), + dropvel_z_(0.1), + radius_(0.1), + interface_thickness_(0.0025) +{ + // check size and populate + if (params.size() != 8 && !params.empty()) + throw std::runtime_error("Realm::setup_initial_conditions: " + "droplet (velocity) requires 8 params: 3 " + "components of droplet position, 3 " + "components of droplet velocity, droplet " + "radius, and interface thickness"); + if (!params.empty()) { + droppos_x_ = params[0]; + droppos_y_ = params[1]; + droppos_z_ = params[2]; + dropvel_x_ = params[3]; + dropvel_y_ = params[4]; + dropvel_z_ = params[5]; + radius_ = params[6]; + interface_thickness_ = params[7]; + } +} + +void +DropletVelocityAuxFunction::do_evaluate( + const double* coords, + const double /*time*/, + const unsigned spatialDimension, + const unsigned numPoints, + double* fieldPtr, + const unsigned fieldSize, + const unsigned /*beginPos*/, + const unsigned /*endPos*/) const +{ + for (unsigned p = 0; p < numPoints; ++p) { + + const double x = coords[0]; + const double y = coords[1]; + const double z = coords[2]; + + auto rad_pos = std::sqrt( + (x - droppos_x_) * (x - droppos_x_) + + (y - droppos_y_) * (y - droppos_y_) + + (z - droppos_z_) * (z - droppos_z_)) - + radius_; + auto vof = -0.5 * (std::erf(rad_pos / interface_thickness_) + 1.0) + 1.0; + + // Approximate average velocity by scaling with vof instead of using density + fieldPtr[0] = vof * dropvel_x_; + fieldPtr[1] = vof * dropvel_y_; + fieldPtr[2] = vof * dropvel_z_; + + fieldPtr += fieldSize; + coords += spatialDimension; + } +} + +} // namespace nalu +} // namespace sierra diff --git a/src/user_functions/SloshingTankVOFAuxFunction.C b/src/user_functions/SloshingTankVOFAuxFunction.C new file mode 100644 index 000000000..58ecea7fc --- /dev/null +++ b/src/user_functions/SloshingTankVOFAuxFunction.C @@ -0,0 +1,71 @@ +// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC +// (NTESS), National Renewable Energy Laboratory, University of Texas Austin, +// Northwest Research Associates. Under the terms of Contract DE-NA0003525 +// with NTESS, the U.S. Government retains certain rights in this software. +// +// This software is released under the BSD 3-clause license. See LICENSE file +// for more details. +// + +#include +#include + +// basic c++ +#include +#include +#include +#include + +namespace sierra { +namespace nalu { + +SloshingTankVOFAuxFunction::SloshingTankVOFAuxFunction( + const std::vector& params) + : AuxFunction(0, 1), + water_level_(0.), + amplitude_(0.1), + kappa_(0.25), + interface_thickness_(0.1) +{ + // check size and populate + if (params.size() != 4 && !params.empty()) + throw std::runtime_error("Realm::setup_initial_conditions: " + "sloshing_tank requires 4 params: water level, " + "amplitude, kappa, and interface thickness"); + if (!params.empty()) { + water_level_ = params[0]; + amplitude_ = params[1]; + kappa_ = params[2]; + interface_thickness_ = params[3]; + } +} + +void +SloshingTankVOFAuxFunction::do_evaluate( + const double* coords, + const double /*time*/, + const unsigned spatialDimension, + const unsigned numPoints, + double* fieldPtr, + const unsigned fieldSize, + const unsigned /*beginPos*/, + const unsigned /*endPos*/) const +{ + for (unsigned p = 0; p < numPoints; ++p) { + + const double x = coords[0]; + const double y = coords[1]; + const double z = coords[2]; + + const double z0 = + water_level_ + amplitude_ * std::exp(-kappa_ * (x * x + y * y)); + fieldPtr[0] = + -0.5 * (std::erf((z - z0) / interface_thickness_) + 1.0) + 1.0; + + fieldPtr += fieldSize; + coords += spatialDimension; + } +} + +} // namespace nalu +} // namespace sierra diff --git a/unit_tests/kernels/UnitTestKernelUtils.h b/unit_tests/kernels/UnitTestKernelUtils.h index 2f134e893..1919ea237 100644 --- a/unit_tests/kernels/UnitTestKernelUtils.h +++ b/unit_tests/kernels/UnitTestKernelUtils.h @@ -1434,13 +1434,15 @@ class VOFKernelHex8Mesh : public TestKernelHex8Mesh dvolumeOfFluidDx_(&meta_->declare_field( stk::topology::NODE_RANK, "volume_of_fluid_gradient")), velocity_( - &meta_->declare_field(stk::topology::NODE_RANK, "velocity")), + &meta_->declare_field(stk::topology::NODE_RANK, "velocity", 2)), density_( &meta_->declare_field(stk::topology::NODE_RANK, "density", 2)), viscosity_( &meta_->declare_field(stk::topology::NODE_RANK, "viscosity")), massFlowRateEdge_(&meta_->declare_field( stk::topology::EDGE_RANK, "mass_flow_rate")), + vofBalancedMassFlowRateEdge_(&meta_->declare_field( + stk::topology::EDGE_RANK, "mass_vof_balanced_flow_rate")), znot_(1.0), amf_(2.0), rhoPrimary_(1000.0), @@ -1463,6 +1465,8 @@ class VOFKernelHex8Mesh : public TestKernelHex8Mesh stk::mesh::put_field_on_mesh(*viscosity_, meta_->universal_part(), nullptr); stk::mesh::put_field_on_mesh( *massFlowRateEdge_, meta_->universal_part(), nullptr); + stk::mesh::put_field_on_mesh( + *vofBalancedMassFlowRateEdge_, meta_->universal_part(), 1, nullptr); } virtual ~VOFKernelHex8Mesh() {} @@ -1491,6 +1495,7 @@ class VOFKernelHex8Mesh : public TestKernelHex8Mesh sierra::nalu::ScalarFieldType* density_{nullptr}; sierra::nalu::ScalarFieldType* viscosity_{nullptr}; sierra::nalu::ScalarFieldType* massFlowRateEdge_{nullptr}; + sierra::nalu::ScalarFieldType* vofBalancedMassFlowRateEdge_{nullptr}; const double znot_; const double amf_; const double rhoPrimary_;