diff --git a/include/NaluParsing.h b/include/NaluParsing.h index ef46f9d6..13cae1c0 100644 --- a/include/NaluParsing.h +++ b/include/NaluParsing.h @@ -189,6 +189,7 @@ struct WallUserData : public UserData { TurbKinEnergy tke_; MixtureFraction mixFrac_; MassFraction massFraction_; + VolumeOfFluid vof_; Emissivity emissivity_; Irradiation irradiation_; Transmissivity transmissivity_; diff --git a/include/SolutionOptions.h b/include/SolutionOptions.h index 38a0ae47..e13f91d2 100644 --- a/include/SolutionOptions.h +++ b/include/SolutionOptions.h @@ -137,7 +137,7 @@ class SolutionOptions double mdotAlgAccumulation_; double mdotAlgInflow_; double mdotAlgOpen_; - + // turbulence model coeffs std::map turbModelConstantMap_; diff --git a/include/VolumeOfFluidEquationSystem.h b/include/VolumeOfFluidEquationSystem.h index 95114e5e..00eaf69f 100644 --- a/include/VolumeOfFluidEquationSystem.h +++ b/include/VolumeOfFluidEquationSystem.h @@ -34,7 +34,11 @@ class VolumeOfFluidEquationSystem : public EquationSystem { VolumeOfFluidEquationSystem( EquationSystems& equationSystems, const bool outputClippingDiag, - const double deltaZClip); + const double deltaVofClip, + const double Fo, + const double cAlpha, + const bool smooth, + const int smoothIter); virtual ~VolumeOfFluidEquationSystem(); void populate_derived_quantities(); @@ -86,15 +90,15 @@ class VolumeOfFluidEquationSystem : public EquationSystem { void sharpen_interface_explicit(); void smooth_vof(); + void smooth_vof_execute(); void compute_interface_normal(); const bool managePNG_; - const bool outputClippingDiag_; - const double deltaVofClip_; ScalarFieldType *vof_; ScalarFieldType *vofSmoothed_; - ScalarFieldType *interfaceNormal_; + ScalarFieldType *smoothedRhs_; + VectorFieldType *interfaceNormal_; VectorFieldType *dvofdx_; ScalarFieldType *vofTmp_; @@ -102,6 +106,17 @@ class VolumeOfFluidEquationSystem : public EquationSystem { ProjectedNodalGradientEquationSystem *projectedNodalGradEqs_; + // clipping + const bool outputClippingDiag_; + const double deltaVofClip_; + + // smoothing and sharepening params + const double Fo_; + const double cAlpha_; + double dxMin_; + const bool smooth_; + const int smoothIter_; + bool isInit_; }; diff --git a/include/kernel/VolumeOfFluidSharpenElemKernel.h b/include/kernel/VolumeOfFluidSharpenElemKernel.h new file mode 100644 index 00000000..7dcbb9df --- /dev/null +++ b/include/kernel/VolumeOfFluidSharpenElemKernel.h @@ -0,0 +1,65 @@ +/*------------------------------------------------------------------------*/ +/* Copyright 2014 Sandia Corporation. */ +/* This software is released under the license detailed */ +/* in the file, LICENSE, which is located in the top-level Nalu */ +/* directory structure */ +/*------------------------------------------------------------------------*/ + +#ifndef VolumeOfFluidSharpenElemKernel_H +#define VolumeOfFluidSharpenElemKernel_H + +#include "kernel/Kernel.h" +#include "FieldTypeDef.h" + +#include +#include + +#include + +namespace sierra { +namespace nalu { + +class TimeIntegrator; +class SolutionOptions; +class MasterElement; +class ElemDataRequests; + +template +class VolumeOfFluidSharpenElemKernel: public Kernel +{ +public: + VolumeOfFluidSharpenElemKernel( + const stk::mesh::BulkData&, + const SolutionOptions&, + ScalarFieldType*, + const double, + ElemDataRequests&); + + virtual ~VolumeOfFluidSharpenElemKernel(); + + /** Execute the kernel within a Kokkos loop and populate the LHS and RHS for + * the linear solve + */ + virtual void execute( + SharedMemView&, + SharedMemView&, + ScratchViews&); + +private: + VolumeOfFluidSharpenElemKernel() = delete; + + VectorFieldType *coordinates_{nullptr}; + VectorFieldType *velocityRTM_{nullptr}; + VectorFieldType *interfaceNormal_{nullptr}; + ScalarFieldType *vofNp1_{nullptr}; + + /// Integration point to node mapping + const double cAlpha_; + const int* ipNodeMap_; + +}; + +} // nalu +} // sierra + +#endif /* VolumeOfFluidSharpenElemKernel_H */ diff --git a/include/user_functions/RayleighTaylorMixFracAuxFunction.h b/include/user_functions/RayleighTaylorMixFracAuxFunction.h index 2aff1fac..f50a6088 100644 --- a/include/user_functions/RayleighTaylorMixFracAuxFunction.h +++ b/include/user_functions/RayleighTaylorMixFracAuxFunction.h @@ -39,6 +39,7 @@ class RayleighTaylorMixFracAuxFunction : public AuxFunction double tX_; double yTr_; double dTr_; + double surf_; const double pi_; }; diff --git a/reg_tests/test_files/vofSphereDrop/vofSphereDrop.i b/reg_tests/test_files/vofSphereDrop/vofSphereDrop.i index c81940ee..1a79bd8b 100644 --- a/reg_tests/test_files/vofSphereDrop/vofSphereDrop.i +++ b/reg_tests/test_files/vofSphereDrop/vofSphereDrop.i @@ -51,6 +51,9 @@ realms: name: myV max_iterations: 1 convergence_tolerance: 1.e-2 + activate_smoothing: yes + fourier_number: 0.25 + compression_constant: 0.05 initial_conditions: - constant: ic_1 @@ -149,7 +152,7 @@ realms: - element_source_terms: momentum: [lumped_momentum_time_derivative, advection_diffusion, buoyancy] continuity: [advection] - volume_of_fluid: [vof, sucv_nso] + volume_of_fluid: [vof, sucv_nso, sharpen] - user_constants: gravity: [0.0, -10.0] @@ -169,6 +172,10 @@ realms: - intersected_element - mesh_velocity - density + - volume_of_fluid + - volume_of_fluid_smoothed + - interface_normal + - dvofdx Time_Integrators: - StandardTimeIntegrator: diff --git a/reg_tests/test_files/vofSphereDrop/vofSphereDrop.norm.gold b/reg_tests/test_files/vofSphereDrop/vofSphereDrop.norm.gold index 24bd328c..981c5ae5 100644 --- a/reg_tests/test_files/vofSphereDrop/vofSphereDrop.norm.gold +++ b/reg_tests/test_files/vofSphereDrop/vofSphereDrop.norm.gold @@ -1,20 +1,20 @@ -0.02432116329066863 1 0.005 -0.001966663692513982 2 0.01 -0.0003478067418379743 3 0.015 -0.000128807125321143 4 0.02 -0.0001186665578123386 5 0.025 -6.159589527807734e-05 6 0.03 -2.85084248670952e-05 7 0.035 -1.804536444197644e-05 8 0.04 -1.696814237987382e-05 9 0.045 -1.82415952817005e-05 10 0.05 -2.013158697076528e-05 11 0.055 -2.227717639712913e-05 12 0.06 -2.459260313863444e-05 13 0.065 -2.704331805762574e-05 14 0.07 -2.961771509211273e-05 15 0.075 -3.230282985030123e-05 16 0.08 -3.50955110186631e-05 17 0.085 -3.796971046972843e-05 18 0.09 -4.093875506528796e-05 19 0.095 -4.399335630296404e-05 20 0.1 +0.01447384532257732 1 0.005 +0.01216518334441016 2 0.01 +0.01073768042783061 3 0.015 +0.01051559281054295 4 0.02 +0.01048860402882239 5 0.025 +0.01041243777094482 6 0.03 +0.01035905603141019 7 0.035 +0.01032539481913046 8 0.04 +0.01029914072916844 9 0.045 +0.0102748391713822 10 0.05 +0.01025087973701757 11 0.055 +0.01022698076916539 12 0.06 +0.01020317723200227 13 0.065 +0.01017951831366472 14 0.07 +0.01015602635420198 15 0.075 +0.01013270865108439 16 0.08 +0.01010957223961167 17 0.085 +0.01008660623117398 18 0.09 +0.01006382646893597 19 0.095 +0.01004123400529326 20 0.1 diff --git a/src/EquationSystems.C b/src/EquationSystems.C index 95398c20..09cfef94 100644 --- a/src/EquationSystems.C +++ b/src/EquationSystems.C @@ -104,7 +104,17 @@ void EquationSystems::load(const YAML::Node & y_node) get_if_present_no_default(y_eqsys, "output_clipping_diagnostic", outputClipDiag); double deltaVofClip = 0.0; get_if_present_no_default(y_eqsys, "clipping_delta", deltaVofClip); - eqSys = new VolumeOfFluidEquationSystem(*this, outputClipDiag, deltaVofClip); + // vof smoothing/sharpening; provide defaults... + double fourierNumber = 0.25; + double cAlpha = 0.05; + bool smooth = true; + int smoothIter = 5; + get_if_present_no_default(y_eqsys, "fourier_number", fourierNumber); + get_if_present_no_default(y_eqsys, "compression_constant", cAlpha); + get_if_present_no_default(y_eqsys, "activate_smoothing", smooth); + get_if_present_no_default(y_eqsys, "smoothing_iterations", smoothIter); + eqSys = new VolumeOfFluidEquationSystem(*this, outputClipDiag, deltaVofClip, + fourierNumber, cAlpha, smooth, smoothIter); } else if ( expect_map(y_system, "LowMachEOM", true) ) { y_eqsys = expect_map(y_system, "LowMachEOM", true); diff --git a/src/NaluParsing.C b/src/NaluParsing.C index 2698f0d5..5df6ef6f 100644 --- a/src/NaluParsing.C +++ b/src/NaluParsing.C @@ -821,6 +821,13 @@ namespace YAML wallData.bcDataSpecifiedMap_["mass_fraction"] = true; wallData.bcDataTypeMap_["mass_fraction"] = sierra::nalu::CONSTANT_UD; } + if (node["volume_of_fluid"]) + { + wallData.vof_ = node["volume_of_fluid"].as< + sierra::nalu::VolumeOfFluid>(); + wallData.bcDataSpecifiedMap_["volume_of_fluid"] = true; + wallData.bcDataTypeMap_["volume_of_fluid"] = sierra::nalu::CONSTANT_UD; + } if (node["emissivity"]) { wallData.emissivity_ = node["emissivity"].as(); diff --git a/src/VolumeOfFluidEquationSystem.C b/src/VolumeOfFluidEquationSystem.C index 8ec2941d..9afcf421 100644 --- a/src/VolumeOfFluidEquationSystem.C +++ b/src/VolumeOfFluidEquationSystem.C @@ -1,4 +1,4 @@ -/*------------------------------------------------------------------------*/ +//*------------------------------------------------------------------------*/ /* Copyright 2014 Sandia Corporation. */ /* This software is released under the license detailed */ /* in the file, LICENSE, which is located in the top-level Nalu */ @@ -41,6 +41,7 @@ #include "AssembleElemSolverAlgorithm.h" #include "kernel/VolumeOfFluidElemKernel.h" #include "kernel/VolumeOfFluidSucvNsoElemKernel.h" +#include "kernel/VolumeOfFluidSharpenElemKernel.h" // user function #include "user_functions/RayleighTaylorMixFracAuxFunction.h" @@ -85,18 +86,28 @@ namespace nalu{ VolumeOfFluidEquationSystem::VolumeOfFluidEquationSystem( EquationSystems& eqSystems, const bool outputClippingDiag, - const double deltaVofClip) + const double deltaVofClip, + const double Fo, + const double cAlpha, + const bool smooth, + const int smoothIter) : EquationSystem(eqSystems, "VolumeOfFluidEQS", "volume_of_fluid"), managePNG_(realm_.get_consistent_mass_matrix_png("volume_of_fluid")), - outputClippingDiag_(outputClippingDiag), - deltaVofClip_(deltaVofClip), vof_(NULL), vofSmoothed_(NULL), + smoothedRhs_(NULL), interfaceNormal_(NULL), dvofdx_(NULL), vofTmp_(NULL), assembleNodalGradAlgDriver_(new AssembleNodalGradAlgorithmDriver(realm_, "volume_of_fluid", "dvofdx")), projectedNodalGradEqs_(NULL), + outputClippingDiag_(outputClippingDiag), + deltaVofClip_(deltaVofClip), + Fo_(Fo), + cAlpha_(cAlpha), + dxMin_(1.0e16), + smooth_(smooth), + smoothIter_(smoothIter), isInit_(true) { // extract solver name and solver object @@ -119,6 +130,14 @@ VolumeOfFluidEquationSystem::VolumeOfFluidEquationSystem( if ( managePNG_ ) { manage_projected_nodal_gradient(eqSystems); } + + // output options + if ( smooth_ ) { + NaluEnv::self().naluOutputP0() << "SCS smoothing_iterations: " << smoothIter_ << std::endl; + NaluEnv::self().naluOutputP0() << "fourier_number: " << Fo_ << std::endl; + } + if ( supp_alg_is_requested("sharpen") ) + NaluEnv::self().naluOutputP0() << "compression_constant: " << cAlpha_ << std::endl; } //-------------------------------------------------------------------------- @@ -156,13 +175,17 @@ VolumeOfFluidEquationSystem::register_nodal_fields( stk::mesh::put_field_on_mesh(*vof_, *part, nullptr); realm_.augment_restart_variable_list("volume_of_fluid"); - // smoothed - vofSmoothed_ = &(meta_data.declare_field(stk::topology::NODE_RANK, "volume_of_fluid_smoothed")); - stk::mesh::put_field_on_mesh(*vofSmoothed_, *part, nullptr); + // smoothed field and smoothed rhs + if ( smooth_ ) { + vofSmoothed_ = &(meta_data.declare_field(stk::topology::NODE_RANK, "volume_of_fluid_smoothed")); + stk::mesh::put_field_on_mesh(*vofSmoothed_, *part, nullptr); + smoothedRhs_ = &(meta_data.declare_field(stk::topology::NODE_RANK, "smoothed_rhs")); + stk::mesh::put_field_on_mesh(*smoothedRhs_, *part, nullptr); + } // normal - interfaceNormal_ = &(meta_data.declare_field(stk::topology::NODE_RANK, "interface_normal")); - stk::mesh::put_field_on_mesh(*interfaceNormal_, *part, nullptr); + interfaceNormal_ = &(meta_data.declare_field(stk::topology::NODE_RANK, "interface_normal")); + stk::mesh::put_field_on_mesh(*interfaceNormal_, *part, nDim, nullptr); // projected nodal gradient dvofdx_ = &(meta_data.declare_field(stk::topology::NODE_RANK, "dvofdx")); @@ -197,7 +220,8 @@ VolumeOfFluidEquationSystem::register_interior_algorithm( // types of algorithms const AlgorithmType algType = INTERIOR; - ScalarFieldType &vofNp1 = vofSmoothed_->field_of_state(stk::mesh::StateNone); + ScalarFieldType &vofNp1 = (smooth_) ? vofSmoothed_->field_of_state(stk::mesh::StateNone) + : vof_->field_of_state(stk::mesh::StateNP1); VectorFieldType &dvofdxNone = dvofdx_->field_of_state(stk::mesh::StateNone); // non-solver; contribution to projected nodal gradient; allow for element-based shifted @@ -245,6 +269,10 @@ VolumeOfFluidEquationSystem::register_interior_algorithm( build_topo_kernel_if_requested (partTopo, *this, activeKernels, "nso", realm_.bulk_data(), *realm_.solutionOptions_, vof_, 0.0, 1.0, dataPreReqs); + + build_topo_kernel_if_requested + (partTopo, *this, activeKernels, "sharpen", + realm_.bulk_data(), *realm_.solutionOptions_, vof_, cAlpha_, dataPreReqs); report_invalid_supp_alg_names(); report_built_supp_alg_names(); @@ -264,7 +292,8 @@ VolumeOfFluidEquationSystem::register_inflow_bc( // algorithm type const AlgorithmType algType = INFLOW; - ScalarFieldType &vofNp1 = vofSmoothed_->field_of_state(stk::mesh::StateNone); + ScalarFieldType &vofNp1 = (smooth_) ? vofSmoothed_->field_of_state(stk::mesh::StateNone) + : vof_->field_of_state(stk::mesh::StateNP1); VectorFieldType &dvofdxNone = dvofdx_->field_of_state(stk::mesh::StateNone); stk::mesh::MetaData &meta_data = realm_.meta_data(); @@ -356,7 +385,8 @@ VolumeOfFluidEquationSystem::register_open_bc( // algorithm type const AlgorithmType algType = OPEN; - ScalarFieldType &vofNp1 = vofSmoothed_->field_of_state(stk::mesh::StateNone); + ScalarFieldType &vofNp1 = (smooth_) ? vofSmoothed_->field_of_state(stk::mesh::StateNone) + : vof_->field_of_state(stk::mesh::StateNP1); VectorFieldType &dvofdxNone = dvofdx_->field_of_state(stk::mesh::StateNone); // non-solver; dvofdx; allow for element-based shifted @@ -388,9 +418,60 @@ VolumeOfFluidEquationSystem::register_wall_bc( // algorithm type const AlgorithmType algType = WALL; - // np1 - ScalarFieldType &vofNp1 = vofSmoothed_->field_of_state(stk::mesh::StateNone); + ScalarFieldType &vofNp1 = (smooth_) ? vofSmoothed_->field_of_state(stk::mesh::StateNone) + : vof_->field_of_state(stk::mesh::StateNP1); VectorFieldType &dvofdxNone = dvofdx_->field_of_state(stk::mesh::StateNone); + + stk::mesh::MetaData & meta_data = realm_.meta_data(); + + // extract the value for [optional] user specified vof + WallUserData userData = wallBCData.userData_; + std::string vofName = "volume_of_fluid"; + if ( bc_data_specified(userData, vofName) ) { + + // FIXME: Generalize for constant vs function + + ScalarFieldType &realVofNp1 = vof_->field_of_state(stk::mesh::StateNP1); + + // register boundary data; vof_bc + ScalarFieldType *theBcField = &(meta_data.declare_field(stk::topology::NODE_RANK, "vof_bc")); + stk::mesh::put_field_on_mesh(*theBcField, *part, nullptr); + + // extract data + std::vector userSpec(1); + VolumeOfFluid vof = userData.vof_; + userSpec[0] = vof.vof_; + + // new it + ConstantAuxFunction *theAuxFunc = new ConstantAuxFunction(0, 1, userSpec); + + // bc data alg + AuxFunctionAlgorithm *auxAlg + = new AuxFunctionAlgorithm(realm_, part, + theBcField, theAuxFunc, + stk::topology::NODE_RANK); + bcDataAlg_.push_back(auxAlg); + + // copy vof_bc to vof np1... + CopyFieldAlgorithm *theCopyAlg + = new CopyFieldAlgorithm(realm_, part, + theBcField, &realVofNp1, + 0, 1, + stk::topology::NODE_RANK); + bcDataMapAlg_.push_back(theCopyAlg); + + // Dirichlet bc + std::map::iterator itd = + solverAlgDriver_->solverDirichAlgMap_.find(algType); + if ( itd == solverAlgDriver_->solverDirichAlgMap_.end() ) { + DirichletBC *theAlg + = new DirichletBC(realm_, this, part, &realVofNp1, theBcField, 0, 1); + solverAlgDriver_->solverDirichAlgMap_[algType] = theAlg; + } + else { + itd->second->partVec_.push_back(part); + } + } // non-solver; dvofdx; allow for element-based shifted if ( !managePNG_ ) { @@ -406,7 +487,7 @@ VolumeOfFluidEquationSystem::register_wall_bc( } } } - + //-------------------------------------------------------------------------- //-------- register_symmetry_bc -------------------------------------------- //-------------------------------------------------------------------------- @@ -420,8 +501,8 @@ VolumeOfFluidEquationSystem::register_symmetry_bc( // algorithm type const AlgorithmType algType = SYMMETRY; - // np1 - ScalarFieldType &vofNp1 = vofSmoothed_->field_of_state(stk::mesh::StateNone); + ScalarFieldType &vofNp1 = (smooth_) ? vofSmoothed_->field_of_state(stk::mesh::StateNone) + : vof_->field_of_state(stk::mesh::StateNP1); VectorFieldType &dvofdxNone = dvofdx_->field_of_state(stk::mesh::StateNone); // non-solver; dvofdx; allow for element-based shifted @@ -450,10 +531,9 @@ VolumeOfFluidEquationSystem::register_overset_bc() UpdateOversetFringeAlgorithmDriver* theAlg = new UpdateOversetFringeAlgorithmDriver(realm_); // Perform fringe updates before all equation system solves equationSystems_.preIterAlgDriver_.push_back(theAlg); - theAlg->fields_.push_back( std::unique_ptr(new OversetFieldData(vof_,1,1))); - + if ( realm_.has_mesh_motion() ) { UpdateOversetFringeAlgorithmDriver* theAlgPost = new UpdateOversetFringeAlgorithmDriver(realm_,false); // Perform fringe updates after all equation system solves (ideally on the post_time_step) @@ -668,20 +748,24 @@ VolumeOfFluidEquationSystem::update_and_clip() void VolumeOfFluidEquationSystem::compute_interface_normal() { - // interface normal is generally compiuted based on vofSmoothed_ - stk::mesh::MetaData & meta_data = realm_.meta_data(); - stk::mesh::BulkData & bulk_data = realm_.bulk_data(); + stk::mesh::MetaData & metaData = realm_.meta_data(); + stk::mesh::BulkData & bulkData = realm_.bulk_data(); - const int nDim = meta_data.spatial_dimension(); + const int nDim = metaData.spatial_dimension(); const double small = 1.0e-16; // extract nodal fields - VectorFieldType *coordinates = meta_data.get_field(stk::topology::NODE_RANK, realm_.get_coordinates_name()); - - ScalarFieldType *dualNodalVolume = meta_data.get_field(stk::topology::NODE_RANK, "dual_nodal_volume"); - + VectorFieldType *coordinates + = metaData.get_field(stk::topology::NODE_RANK, realm_.get_coordinates_name()); + ScalarFieldType *dualNodalVolume + = metaData.get_field(stk::topology::NODE_RANK, "dual_nodal_volume"); + + // interface normal is generally computed based on vofSmoothed_; user defines this + ScalarFieldType &vofNp1 = (smooth_) ? vofSmoothed_->field_of_state(stk::mesh::StateNone) + : vof_->field_of_state(stk::mesh::StateNP1); + // zero assembled normal - field_fill( meta_data, bulk_data, 0.0, *interfaceNormal_, realm_.get_activate_aura()); + field_fill( metaData, bulkData, 0.0, *interfaceNormal_, realm_.get_activate_aura()); // integration point data that is fixed std::vector dvofdxIp(nDim); @@ -689,16 +773,16 @@ VolumeOfFluidEquationSystem::compute_interface_normal() // fields to gather std::vector ws_coordinates; - std::vector ws_vofSmoothed; - std::vector ws_dual_volume; - std::vector ws_scv_volume; + std::vector ws_vofNp1; + std::vector ws_dualVolume; + std::vector ws_scVolume; std::vector ws_dndx; std::vector ws_deriv; std::vector ws_det_j; // select locally owned where vof is defined; exclude inactive block - stk::mesh::Selector s_locally_owned_union = meta_data.locally_owned_part() - & stk::mesh::selectField(*vofSmoothed_) + stk::mesh::Selector s_locally_owned_union = metaData.locally_owned_part() + & stk::mesh::selectField(vofNp1) & !(realm_.get_inactive_selector()); stk::mesh::BucketVector const& elem_buckets = @@ -718,18 +802,18 @@ VolumeOfFluidEquationSystem::compute_interface_normal() // algorithm related ws_coordinates.resize(nodesPerElement*nDim); - ws_vofSmoothed.resize(nodesPerElement); - ws_dual_volume.resize(nodesPerElement); - ws_scv_volume.resize(numScvIp); + ws_vofNp1.resize(nodesPerElement); + ws_dualVolume.resize(nodesPerElement); + ws_scVolume.resize(numScvIp); ws_dndx.resize(nDim*numScvIp*nodesPerElement); ws_deriv.resize(nDim*numScvIp*nodesPerElement); ws_det_j.resize(numScvIp); // pointers double *p_coordinates = &ws_coordinates[0]; - double *p_vofSmoothed = &ws_vofSmoothed[0]; - double *p_dual_volume = &ws_dual_volume[0]; - double *p_scv_volume = &ws_scv_volume[0]; + double *p_vofNp1 = &ws_vofNp1[0]; + double *p_dualVolume = &ws_dualVolume[0]; + double *p_scVolume = &ws_scVolume[0]; double *p_dndx = &ws_dndx[0]; for ( stk::mesh::Bucket::size_type k = 0 ; k < length ; ++k ) { @@ -750,8 +834,8 @@ VolumeOfFluidEquationSystem::compute_interface_normal() const double * coords = stk::mesh::field_data(*coordinates, node); // gather scalars - p_vofSmoothed[ni] = *stk::mesh::field_data(*vofSmoothed_, node); - p_dual_volume[ni] = *stk::mesh::field_data(*dualNodalVolume, node); + p_vofNp1[ni] = *stk::mesh::field_data(vofNp1, node); + p_dualVolume[ni] = *stk::mesh::field_data(*dualNodalVolume, node); // gather vectors const int offSet = ni*nDim; @@ -761,11 +845,11 @@ VolumeOfFluidEquationSystem::compute_interface_normal() } // compute geometry - double scv_error = 0.0; - meSCV->determinant(1, &p_coordinates[0], &p_scv_volume[0], &scv_error); + double scvError = 0.0; + meSCV->determinant(1, &p_coordinates[0], &p_scVolume[0], &scvError); // compute dndx - meSCV->grad_op(1, &p_coordinates[0], &p_dndx[0], &ws_deriv[0], &ws_det_j[0], &scv_error); + meSCV->grad_op(1, &p_coordinates[0], &p_dndx[0], &ws_deriv[0], &ws_det_j[0], &scvError); for ( int ip = 0; ip < numScvIp; ++ip ) { @@ -777,9 +861,9 @@ VolumeOfFluidEquationSystem::compute_interface_normal() // compute gradient for ( int ic = 0; ic < nodesPerElement; ++ic ) { const int offSetDnDx = nDim*nodesPerElement*ip + ic*nDim; - const double nodalVofSmoothed = p_vofSmoothed[ic]; + const double nodalVof = p_vofNp1[ic]; for ( int j = 0; j < nDim; ++j ) { - p_dvofdxIp[j] += p_dndx[offSetDnDx+j]*nodalVofSmoothed; + p_dvofdxIp[j] += p_dndx[offSetDnDx+j]*nodalVof; } } @@ -797,7 +881,7 @@ VolumeOfFluidEquationSystem::compute_interface_normal() // pointers to real data double * interfaceNormal = stk::mesh::field_data(*interfaceNormal_, node); - const double volumeFac = p_scv_volume[ip]/p_dual_volume[nn]; + const double volumeFac = p_scVolume[ip]/p_dualVolume[nn]; for ( int j = 0; j < nDim; ++j ) interfaceNormal[j] += p_dvofdxIp[j]/dvofMag*volumeFac; } @@ -805,27 +889,204 @@ VolumeOfFluidEquationSystem::compute_interface_normal() } // parallel sum - stk::mesh::parallel_sum(bulk_data, {interfaceNormal_}); + stk::mesh::parallel_sum(bulkData, {interfaceNormal_}); // periodic assemble if ( realm_.hasPeriodic_) { - const unsigned fieldSize = 1; - realm_.periodic_field_update(interfaceNormal_, fieldSize); + realm_.periodic_field_update(interfaceNormal_, nDim); } + // overset update + if ( realm_.hasOverset_ ) { + realm_.overset_orphan_node_field_update(interfaceNormal_, 1, nDim); + } } + //-------------------------------------------------------------------------- //-------- smooth_vof ------------------------------------------------------ //-------------------------------------------------------------------------- void VolumeOfFluidEquationSystem::smooth_vof() { - // for now, just copy vof to vofSmoothed - NaluEnv::self().naluOutputP0() << "VOF::smooth_interface() TBD" << std::endl; - ScalarFieldType &vofNp1 = vof_->field_of_state(stk::mesh::StateN); - ScalarFieldType &vofSmoothed = vofSmoothed_->field_of_state(stk::mesh::StateNone); - field_copy(realm_.meta_data(), realm_.bulk_data(), vofNp1, vofSmoothed, realm_.get_activate_aura()); + + // may not want to smooth + if ( !smooth_) + return; + + // otherwise, proceed + stk::mesh::MetaData & metaData = realm_.meta_data(); + stk::mesh::BulkData & bulkData = realm_.bulk_data(); + + // Copy vof to vofSmoothed + field_copy(metaData, bulkData, *vof_, *vofSmoothed_, realm_.get_activate_aura()); + + // fixed set of iterations... + for( int k = 0; k < smoothIter_; ++k) { + smooth_vof_execute(); + // vofSmoothed = Fo_*dxMin_*dxMin_*smoothedRhs_ + 1.0*vofSmoothed_ + field_axpby(metaData, bulkData, Fo_*dxMin_*dxMin_, *smoothedRhs_, 1.0, *vofSmoothed_, realm_.get_activate_aura()); + } + +} + +//-------------------------------------------------------------------------- +//-------- smooth_vof_execute ---------------------------------------------- +//-------------------------------------------------------------------------- +void +VolumeOfFluidEquationSystem::smooth_vof_execute() +{ + // compute smoothedRhs = (Fo*dx*dx) d^2(vof)/dxj^2 using a low-order lumped projection + // leave off smoothing factors until axpby and fold in the dxMin_ here.. + double l_dxMin = 1.0e16; + + stk::mesh::MetaData & metaData = realm_.meta_data(); + stk::mesh::BulkData & bulkData = realm_.bulk_data(); + + const int nDim = metaData.spatial_dimension(); + + // extract nodal fields + VectorFieldType *coordinates + = metaData.get_field(stk::topology::NODE_RANK, realm_.get_coordinates_name()); + ScalarFieldType *dualNodalVolume + = metaData.get_field(stk::topology::NODE_RANK, "dual_nodal_volume"); + + // zero it + field_fill(metaData, bulkData, 0.0, *smoothedRhs_, realm_.get_activate_aura()); + + // geometry related to populate + std::vector ws_vof; + std::vector ws_dualVolume; + std::vector ws_coordinates; + + // define some common selectors + stk::mesh::Selector s_locally_owned_union = metaData.locally_owned_part() + & stk::mesh::selectField(*vofSmoothed_) + & !(realm_.get_inactive_selector()); + + stk::mesh::BucketVector const& elem_buckets = + realm_.get_buckets( stk::topology::ELEMENT_RANK, s_locally_owned_union ); + + std::vector ws_scsAreav; + std::vector ws_dndx; + std::vector ws_deriv; + std::vector ws_detJ; + + for ( const stk::mesh::Bucket* bucket_ptr : elem_buckets ) { + const stk::mesh::Bucket & b = *bucket_ptr ; + const stk::mesh::Bucket::size_type length = b.size(); + + // extract master element + MasterElement *meSCS = sierra::nalu::MasterElementRepo::get_surface_master_element(b.topology()); + + // extract master element specifics + const int nodesPerElement = meSCS->nodesPerElement_; + const int numScsIp = meSCS->numIntPoints_; + const int *lrscv = meSCS->adjacentNodes(); + + // algorithm related + ws_coordinates.resize(nodesPerElement*nDim); + ws_vof.resize(nodesPerElement); + ws_dualVolume.resize(nodesPerElement); + ws_scsAreav.resize(numScsIp*nDim); + ws_dndx.resize(nDim*numScsIp*nodesPerElement); + ws_deriv.resize(nDim*numScsIp*nodesPerElement); + ws_detJ.resize(numScsIp); + + // pointers + double *p_coordinates = &ws_coordinates[0]; + double *p_vof = &ws_vof[0]; + double *p_dualVolume = &ws_dualVolume[0]; + double *p_scsAreav = &ws_scsAreav[0]; + double *p_dndx = &ws_dndx[0]; + + for ( stk::mesh::Bucket::size_type k = 0 ; k < length ; ++k ) { + + //=============================================== + // gather nodal data; this is how we do it now.. + //=============================================== + stk::mesh::Entity const * node_rels = b.begin_nodes(k); + int num_nodes = b.num_nodes(k); + + // sanity check on num nodes + ThrowAssert( num_nodes == nodesPerElement ); + + for ( int ni = 0; ni < num_nodes; ++ni ) { + stk::mesh::Entity node = node_rels[ni]; + + // pointers to real data + const double * coords = stk::mesh::field_data(*coordinates, node ); + + // gather scalars + p_vof[ni] = *stk::mesh::field_data(*vofSmoothed_, node ); + p_dualVolume[ni] = *stk::mesh::field_data(*dualNodalVolume, node ); + + // gather vectors + const int niNdim = ni*nDim; + for ( int j=0; j < nDim; ++j ) { + p_coordinates[niNdim+j] = coords[j]; + } + } + + // compute geometry + double scsError = 0.0; + meSCS->determinant(1, &p_coordinates[0], &p_scsAreav[0], &scsError); + meSCS->grad_op(1, &p_coordinates[0], &p_dndx[0], &ws_deriv[0], &ws_detJ[0], &scsError); + + for ( int ip = 0; ip < numScsIp; ++ip ) { + + // left and right nodes for this ip + const int il = lrscv[2*ip]; + const int ir = lrscv[2*ip+1]; + + // determine dx; edge distance magnitude + double dx = 0.0; + for ( int j = 0; j < nDim; ++j ) { + const double dxj = p_coordinates[ir*nDim+j] - p_coordinates[il*nDim+j]; + dx += dxj*dxj; + } + l_dxMin = std::min(l_dxMin, std::sqrt(dx)); + + stk::mesh::Entity nodeL = node_rels[il]; + stk::mesh::Entity nodeR = node_rels[ir]; + + // pointer to fields to assemble + double *smoothedRhsL = stk::mesh::field_data(*smoothedRhs_, nodeL ); + double *smoothedRhsR = stk::mesh::field_data(*smoothedRhs_, nodeR ); + + double qDiff = 0.0; + for ( int ic = 0; ic < nodesPerElement; ++ic ) { + + double lhsfacDiff = 0.0; + const int offSetDnDx = nDim*nodesPerElement*ip + ic*nDim; + for ( int j = 0; j < nDim; ++j ) { + lhsfacDiff += -p_dndx[offSetDnDx+j]*p_scsAreav[ip*nDim+j]; + } + qDiff += lhsfacDiff*p_vof[ic]; + } + + *smoothedRhsL -= qDiff/ws_dualVolume[il]; + *smoothedRhsR += qDiff/ws_dualVolume[ir]; + } + } + } + + // parallel sum + stk::mesh::parallel_sum(bulkData, {smoothedRhs_}); + + // periodic update + if ( realm_.hasPeriodic_) { + realm_.periodic_field_update(smoothedRhs_, 1); + } + + // overset update + if ( realm_.hasOverset_ ) { + realm_.overset_orphan_node_field_update(smoothedRhs_, 1, 1); + } + + // compute min + stk::ParallelMachine comm = NaluEnv::self().parallel_comm(); + stk::all_reduce_min(comm, &l_dxMin, &dxMin_, 1); } //-------------------------------------------------------------------------- @@ -834,7 +1095,7 @@ VolumeOfFluidEquationSystem::smooth_vof() void VolumeOfFluidEquationSystem::sharpen_interface_explicit() { - NaluEnv::self().naluOutputP0() << "VOF::sharpen_interface() currently in situ via the linear system" << std::endl; + NaluEnv::self().naluOutputP0() << "VOF::sharpen_interface_explicit() inative" << std::endl; } //-------------------------------------------------------------------------- diff --git a/src/kernel/VolumeOfFluidSharpenElemKernel.C b/src/kernel/VolumeOfFluidSharpenElemKernel.C new file mode 100644 index 00000000..103387b0 --- /dev/null +++ b/src/kernel/VolumeOfFluidSharpenElemKernel.C @@ -0,0 +1,122 @@ +/*------------------------------------------------------------------------*/ +/* Copyright 2014 Sandia Corporation. */ +/* This software is released under the license detailed */ +/* in the file, LICENSE, which is located in the top-level Nalu */ +/* directory structure */ +/*------------------------------------------------------------------------*/ + +#include "kernel/VolumeOfFluidSharpenElemKernel.h" +#include "AlgTraits.h" +#include "master_element/MasterElement.h" +#include "TimeIntegrator.h" +#include "SolutionOptions.h" + +// template and scratch space +#include "BuildTemplates.h" +#include "ScratchViews.h" + +// stk_mesh/base/fem +#include +#include +#include +#include + +namespace sierra { +namespace nalu { + +template +VolumeOfFluidSharpenElemKernel::VolumeOfFluidSharpenElemKernel( + const stk::mesh::BulkData& bulkData, + const SolutionOptions& solnOpts, + ScalarFieldType* vof, + const double cAlpha, + ElemDataRequests& dataPreReqs) + : Kernel(), + cAlpha_(cAlpha), + ipNodeMap_(sierra::nalu::MasterElementRepo::get_volume_master_element(AlgTraits::topo_)->ipNodeMap()) +{ + // save off fields + const stk::mesh::MetaData& metaData = bulkData.mesh_meta_data(); + + vofNp1_ = &(vof->field_of_state(stk::mesh::StateNP1)); + interfaceNormal_ = metaData.get_field( + stk::topology::NODE_RANK, "interface_normal"); + std::string velocity_name = solnOpts.does_mesh_move() ? "velocity_rtm" : "velocity"; + velocityRTM_ = metaData.get_field( + stk::topology::NODE_RANK, velocity_name); + coordinates_ = metaData.get_field( + stk::topology::NODE_RANK, solnOpts.get_coordinates_name()); + + MasterElement *meSCV = sierra::nalu::MasterElementRepo::get_volume_master_element(AlgTraits::topo_); + + // add master elements + dataPreReqs.add_cvfem_volume_me(meSCV); + + // fields and data + dataPreReqs.add_coordinates_field(*coordinates_, AlgTraits::nDim_, CURRENT_COORDINATES); + dataPreReqs.add_gathered_nodal_field(*velocityRTM_, AlgTraits::nDim_); + dataPreReqs.add_gathered_nodal_field(*interfaceNormal_, AlgTraits::nDim_); + dataPreReqs.add_gathered_nodal_field(*vofNp1_, 1); + dataPreReqs.add_master_element_call(SCV_VOLUME, CURRENT_COORDINATES); + dataPreReqs.add_master_element_call(SCV_GRAD_OP, CURRENT_COORDINATES); +} + +template +VolumeOfFluidSharpenElemKernel::~VolumeOfFluidSharpenElemKernel() +{} + +template +void +VolumeOfFluidSharpenElemKernel::execute( + SharedMemView& lhs, + SharedMemView&rhs, + ScratchViews& scratchViews) +{ + NALU_ALIGNED DoubleType w_magU[AlgTraits::nodesPerElement_]; + + SharedMemView& v_vofNp1 = scratchViews.get_scratch_view_1D(*vofNp1_); + SharedMemView& v_vrtm = scratchViews.get_scratch_view_2D(*velocityRTM_); + SharedMemView& v_interfaceNormal = scratchViews.get_scratch_view_2D(*interfaceNormal_); + SharedMemView& v_scVolume = scratchViews.get_me_views(CURRENT_COORDINATES).scv_volume; + SharedMemView& v_dndx = scratchViews.get_me_views(CURRENT_COORDINATES).dndx_scv; + + // determine nodal velocity magnitude + for ( int ic = 0; ic < AlgTraits::nodesPerElement_; ++ic ) { + DoubleType magU = 0.0; + for ( int k = 0; k < AlgTraits::nDim_; ++k ) { + magU += v_vrtm(ic,k)*v_vrtm(ic,k); + } + w_magU[ic] = stk::math::sqrt(magU); + } + + for ( int ip = 0; ip < AlgTraits::numScvIp_; ++ip ) { + + // nearest node to ip + const int nearestNode = ipNodeMap_[ip]; + const DoubleType scV = v_scVolume(ip); + + DoubleType sharpen = 0.0; + for ( int ic = 0; ic < AlgTraits::nodesPerElement_; ++ic ) { + DoubleType sum = 0.0; + DoubleType lhsOne = 0.0; + DoubleType lhsTwo = 0.0; + const DoubleType vofIc = v_vofNp1(ic); + for ( int j = 0; j < AlgTraits::nDim_; ++j ) { + const DoubleType ucj = cAlpha_*w_magU[ic]*v_interfaceNormal(ic,j); + lhsOne += ucj*v_dndx(ip,ic,j); + lhsTwo -= 2.0*ucj*vofIc*v_dndx(ip,ic,j); + sum += ucj*v_dndx(ip,ic,j); + } + sharpen += sum*(vofIc*(1.0-vofIc)); + lhs(nearestNode,ic) += (lhsOne*1.0+lhsTwo*0.0)*scV; + } + + // assemble rhs + rhs(nearestNode) -= sharpen*scV; + } +} + +INSTANTIATE_KERNEL(VolumeOfFluidSharpenElemKernel); + +} // nalu +} // sierra diff --git a/src/user_functions/RayleighTaylorMixFracAuxFunction.C b/src/user_functions/RayleighTaylorMixFracAuxFunction.C index 6b4aae34..2e00717d 100644 --- a/src/user_functions/RayleighTaylorMixFracAuxFunction.C +++ b/src/user_functions/RayleighTaylorMixFracAuxFunction.C @@ -24,6 +24,7 @@ RayleighTaylorMixFracAuxFunction::RayleighTaylorMixFracAuxFunction( tX_(1.0), yTr_(1.0), dTr_(0.20), + surf_(1.0), pi_(acos(-1.0)) { // extract the params - if they are supplied (optional) @@ -35,6 +36,8 @@ RayleighTaylorMixFracAuxFunction::RayleighTaylorMixFracAuxFunction( yTr_ = theParams[2]; if ( theParams.size() > 3 ) dTr_ = theParams[3]; + if ( theParams.size() > 4 ) + surf_ = theParams[4]; } void @@ -67,7 +70,7 @@ RayleighTaylorMixFracAuxFunction::do_evaluate( value = 1.0; } else { - value = 1.0/2.0*(1.0 - sin(pi_*yy/dTr_)); + value = surf_*1.0/2.0*(1.0 - sin(pi_*yy/dTr_)); } fieldPtr[0] = value;