diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index eca91577..1f402fde 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -27,6 +27,7 @@ set(SOURCES_WPENALTY set(SOURCES ${CMAKE_SOURCE_DIR}/source/KKStools.cc ${CMAKE_SOURCE_DIR}/source/main.cc + ${CMAKE_SOURCE_DIR}/source/QuatFaceCoeffs.cc ${CMAKE_SOURCE_DIR}/source/AdaptMovingFrame.cc ${CMAKE_SOURCE_DIR}/source/computeQDiffs.cc ${CMAKE_SOURCE_DIR}/source/ThreeArgsInterpolationType.cc diff --git a/source/QuatFACOps.cc b/source/QuatFACOps.cc index 2f8c951e..a266f037 100644 --- a/source/QuatFACOps.cc +++ b/source/QuatFACOps.cc @@ -2,7 +2,6 @@ // UT-Battelle, LLC. // Produced at the Lawrence Livermore National Laboratory and // the Oak Ridge National Laboratory -// Written by M.R. Dorr, J.-L. Fattebert and M.E. Wickett // LLNL-CODE-747500 // All rights reserved. // This file is part of AMPE. @@ -681,58 +680,6 @@ void QuatFACOps::initializeOperatorState( } } - -void QuatFACOps::computeFaceCoefs(const double epsilon_q, - const int diffusion_coef_id, - const int grad_q_id, - const double gradient_floor, - const std::string grad_floor_type, - const int face_coef_id) // output -{ - - // Check for negative diffusion coefficients. We don't like them. - // Zero is ok since epsilon^2 is added later - - math::HierarchySideDataOpsReal sideops(d_hierarchy, d_ln_min, - d_ln_max - 1); - - double diffusion_coef_min = sideops.min(diffusion_coef_id); - - if ((diffusion_coef_min + epsilon_q * epsilon_q) < 0.) { - TBOX_ERROR(d_object_name << ": Negative diffusion coefficient passed to " - "computeFaceCoefs()."); - } - - for (int ln = d_ln_min; ln <= d_ln_max; ++ln) { - std::shared_ptr level = d_hierarchy->getPatchLevel(ln); - - for (hier::PatchLevel::Iterator pi(level->begin()); pi != level->end(); - pi++) { - std::shared_ptr patch = *pi; - - std::shared_ptr > diffusion_coef_data( - SAMRAI_SHARED_PTR_CAST, hier::PatchData>( - patch->getPatchData(diffusion_coef_id))); - std::shared_ptr > grad_q_data( - SAMRAI_SHARED_PTR_CAST, hier::PatchData>( - patch->getPatchData(grad_q_id))); - std::shared_ptr > face_coef_data( - SAMRAI_SHARED_PTR_CAST, hier::PatchData>( - patch->getPatchData(face_coef_id))); - - assert(diffusion_coef_data->getDepth() == 1); - computeFaceCoefsOnPatch(*patch, epsilon_q, *diffusion_coef_data, - *grad_q_data, *face_coef_data, gradient_floor, - grad_floor_type); - } - } -#if 0 - double norm = sideops.L1Norm( face_coef_id ); - tbox::plog<<"L1 Norm face_coef_id="<& data) } -void QuatFACOps::setOperatorCoefficients( - const double gamma, const double epsilon_q, const int mobility_id, - const int mobility_deriv_id, const int diffusion_coef_id, - const int face_coef_deriv_id, const int grad_q_id, const int q_id, - const double gradient_floor, const std::string grad_floor_type) +void QuatFACOps::setOperatorCoefficients(const double gamma, + const int mobility_id, + const int mobility_deriv_id, + const int face_coef_id, + const int face_coef_deriv_id, + const int grad_q_id, const int q_id) { - assert(epsilon_q > 0.); - d_gamma = gamma; // Check for non-positive mobility @@ -801,11 +747,6 @@ void QuatFACOps::setOperatorCoefficients( "setOperatorCoefficients()."); } - // Compute the face coefficients - - computeFaceCoefs(epsilon_q, diffusion_coef_id, grad_q_id, gradient_floor, - grad_floor_type, d_face_coef_id); - for (int ln = d_ln_min; ln <= d_ln_max; ++ln) { std::shared_ptr level = d_hierarchy->getPatchLevel(ln); @@ -868,7 +809,7 @@ void QuatFACOps::setOperatorCoefficients( // Set the matrix coefficients d_quat_level_solver[ln - d_ln_min]->setMatrixCoefficients(d_gamma, d_sqrt_m_id, - d_face_coef_id); + face_coef_id); } } @@ -1373,81 +1314,7 @@ void QuatFACOps::checkFluxPatchDataIndex() const } } - -/* -******************************************************************* -* * -* AMR-unaware patch-centered computational kernels. * -* * -******************************************************************* -*/ -void QuatFACOps::computeFaceCoefsOnPatch( - const hier::Patch& patch, const double epsilon_q, - pdat::SideData& diffusion_coef_data, - pdat::SideData& grad_q_data, - pdat::SideData& face_coef_data, // output - const double gradient_floor, const std::string grad_floor_type) const -{ -#ifdef DEBUG_CHECK_ASSERTIONS - assert(patch.inHierarchy()); - assert(diffusion_coef_data.getDepth() == 1); - assert(grad_q_data.getDepth() == NDIM * d_qlen); - assert(face_coef_data.getDepth() == d_qlen); -#endif - - const hier::Box& box = patch.getBox(); - const hier::Index& lower = box.lower(); - const hier::Index& upper = box.upper(); - - const hier::Box& dc_gbox = diffusion_coef_data.getGhostBox(); - const hier::Index& dcglower = dc_gbox.lower(); - const hier::Index& dcgupper = dc_gbox.upper(); - - const hier::Box& gq_gbox = grad_q_data.getGhostBox(); - const hier::Index& gqlower = gq_gbox.lower(); - const hier::Index& gqupper = gq_gbox.upper(); - - const hier::Box& d_gbox = face_coef_data.getGhostBox(); - const hier::Index& dlower = d_gbox.lower(); - const hier::Index& dupper = d_gbox.upper(); - -#if NDIM == 2 - COMPUTE_FACE_COEF2D(lower[0], upper[0], lower[1], upper[1], d_qlen, - epsilon_q, diffusion_coef_data.getPointer(0), - dcglower[0], dcgupper[0] + 1, dcglower[1], dcgupper[1], - diffusion_coef_data.getPointer(1), dcglower[0], - dcgupper[0], dcglower[1], dcgupper[1] + 1, - grad_q_data.getPointer(0), gqlower[0], gqupper[0] + 1, - gqlower[1], gqupper[1], grad_q_data.getPointer(1), - gqlower[0], gqupper[0], gqlower[1], gqupper[1] + 1, - face_coef_data.getPointer(0), dlower[0], dupper[0] + 1, - dlower[1], dupper[1], // output - face_coef_data.getPointer(1), dlower[0], dupper[0], - dlower[1], dupper[1] + 1, // output - gradient_floor, grad_floor_type.c_str()); -#endif -#if NDIM == 3 - COMPUTE_FACE_COEF3D( - lower[0], upper[0], lower[1], upper[1], lower[2], upper[2], d_qlen, - epsilon_q, diffusion_coef_data.getPointer(0), dcglower[0], - dcgupper[0] + 1, dcglower[1], dcgupper[1], dcglower[2], dcgupper[2], - diffusion_coef_data.getPointer(1), dcglower[0], dcgupper[0], dcglower[1], - dcgupper[1] + 1, dcglower[2], dcgupper[2], - diffusion_coef_data.getPointer(2), dcglower[0], dcgupper[0], dcglower[1], - dcgupper[1], dcglower[2], dcgupper[2] + 1, grad_q_data.getPointer(0), - gqlower[0], gqupper[0] + 1, gqlower[1], gqupper[1], gqlower[2], - gqupper[2], grad_q_data.getPointer(1), gqlower[0], gqupper[0], - gqlower[1], gqupper[1] + 1, gqlower[2], gqupper[2], - grad_q_data.getPointer(2), gqlower[0], gqupper[0], gqlower[1], - gqupper[1], gqlower[2], gqupper[2] + 1, face_coef_data.getPointer(0), - dlower[0], dupper[0] + 1, dlower[1], dupper[1], dlower[2], dupper[2], - face_coef_data.getPointer(1), dlower[0], dupper[0], dlower[1], - dupper[1] + 1, dlower[2], dupper[2], face_coef_data.getPointer(2), - dlower[0], dupper[0], dlower[1], dupper[1], dlower[2], dupper[2] + 1, - gradient_floor, grad_floor_type.c_str()); -#endif -} - +// compute face_coef_data using dprime_data and phi_data void QuatFACOps::computeDQuatDPhiFaceCoefsOnPatch( const hier::Patch& patch, pdat::SideData& dprime_data, pdat::CellData& phi_data, @@ -1998,34 +1865,29 @@ void QuatFACOps::evaluateRHS( const int grad_q_copy_id, // for computation of diffusion coefficient const double gradient_floor, const std::string gradient_floor_type, const int mobility_id, const int rotation_index_id, const int q_id, - int rhs_id, const bool use_gradq_for_flux) + int rhs_id) { t_compute_rhs->start(); - if (use_gradq_for_flux) assert(grad_q_id >= 0); + assert(grad_q_id >= 0); assert(grad_q_copy_id >= 0); d_rotation_index_id = rotation_index_id; - const int gq_id = use_gradq_for_flux ? grad_q_id : -1; - // Initialize the output array d_hopscell->setToScalar(rhs_id, 0., false); - computeFaceCoefs(epsilon_q, diffusion_coef_id, grad_q_copy_id, - gradient_floor, gradient_floor_type, - d_face_coef_scratch_id); - for (int ln = d_ln_max; ln >= d_ln_min; ln--) { - accumulateOperatorOnLevel(mobility_id, d_face_coef_scratch_id, q_id, - gq_id, rhs_id, ln, true, false); + accumulateOperatorOnLevel(mobility_id, diffusion_coef_id, q_id, grad_q_id, + rhs_id, ln, true, false); } t_compute_rhs->stop(); } -void QuatFACOps::multiplyDQuatDPhiBlock(const int phase_id, const int out_id) +void QuatFACOps::multiplyDQuatDPhiBlock(const int phase_id, const int out_id, + const int face_coef_id) { // Initialize the output array d_hopscell->setToScalar(out_id, 0., false); @@ -2034,7 +1896,7 @@ void QuatFACOps::multiplyDQuatDPhiBlock(const int phase_id, const int out_id) // mobility derivative times the input phi for (int ln = d_ln_max; ln >= d_ln_min; ln--) { - accumulateOperatorOnLevel(d_m_deriv_id, d_face_coef_id, d_q_local_id, -1, + accumulateOperatorOnLevel(d_m_deriv_id, face_coef_id, d_q_local_id, -1, out_id, ln, false, false); std::shared_ptr level = d_hierarchy->getPatchLevel(ln); diff --git a/source/QuatFACOps.h b/source/QuatFACOps.h index 6fff1ec4..8010b805 100644 --- a/source/QuatFACOps.h +++ b/source/QuatFACOps.h @@ -2,36 +2,11 @@ // UT-Battelle, LLC. // Produced at the Lawrence Livermore National Laboratory and // the Oak Ridge National Laboratory -// Written by M.R. Dorr, J.-L. Fattebert and M.E. Wickett // LLNL-CODE-747500 // All rights reserved. // This file is part of AMPE. // For details, see https://github.com/LLNL/AMPE // Please also read AMPE/LICENSE. -// 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 disclaimer below. -// - Redistributions in binary form must reproduce the above copyright notice, -// this list of conditions and the disclaimer (as noted below) in the -// documentation and/or other materials provided with the distribution. -// - Neither the name of the LLNS/LLNL 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 LAWRENCE LIVERMORE NATIONAL SECURITY, -// LLC, UT BATTELLE, LLC, -// THE U.S. DEPARTMENT OF ENERGY 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. // #ifndef included_QuatFACOps_h #define included_QuatFACOps_h @@ -312,11 +287,11 @@ class QuatFACOps : public SAMRAI::solv::FACOperatorStrategy /* * Set the operator coefficients. */ - void setOperatorCoefficients( - const double time_step, const double epsilon_q, const int mobility_id, - const int mobility_deriv_id, const int diff_coef_id, - const int diff_coef_deriv_id, const int grad_q_id, const int q_id, - const double gradient_floor, const std::string grad_floor_type); + void setOperatorCoefficients(const double time_step, const int mobility_id, + const int mobility_deriv_id, + const int diff_coef_id, + const int diff_coef_deriv_id, + const int grad_q_id, const int q_id); // FACOperatorStrategy virtuals @@ -371,14 +346,15 @@ class QuatFACOps : public SAMRAI::solv::FACOperatorStrategy const double gradient_floor, const std::string gradient_floor_type, const int mobility_id, const int rotations_id, - const int q_id, int rhs_id, const bool use_gradq_for_flux); + const int q_id, int rhs_id); void accumulateOperatorOnLevel(const int mobility_id, const int face_coef_id, const int q_id, const int grad_q_id, int operator_q_id, int ln, bool project, bool error_equation_indicator); - void multiplyDQuatDPhiBlock(const int q_id, const int operator_q_id); + void multiplyDQuatDPhiBlock(const int q_id, const int operator_q_id, + const int face_coeff_id); void applyProjectionOnLevel(const int q_id, const int corr_id, const int err_id, const int ln); diff --git a/source/QuatFaceCoeffs.cc b/source/QuatFaceCoeffs.cc new file mode 100644 index 00000000..e56b494a --- /dev/null +++ b/source/QuatFaceCoeffs.cc @@ -0,0 +1,128 @@ +// Copyright (c) 2018, Lawrence Livermore National Security, LLC and +// UT-Battelle, LLC. +// Produced at the Lawrence Livermore National Laboratory and +// the Oak Ridge National Laboratory +// LLNL-CODE-747500 +// All rights reserved. +// This file is part of AMPE. +// For details, see https://github.com/LLNL/AMPE +// Please also read AMPE/LICENSE. +#include "QuatFaceCoeffs.h" +#include "QuatFort.h" + +#include "SAMRAI/math/HierarchySideDataOpsReal.h" + +using namespace SAMRAI; + +QuatFaceCoeffs::QuatFaceCoeffs(const int qlen, const double epsilon_q, + const double gradient_floor, + const std::string grad_floor_type) + : d_qlen(qlen), + d_epsilon_q(epsilon_q), + d_gradient_floor(gradient_floor), + d_grad_floor_type(grad_floor_type) +{ +} + +// Evaluate coefficient eps_q^2+D_q(phi)/|nabla q| +void QuatFaceCoeffs::computeCoeffs( + const std::shared_ptr hierarchy, + const int diffusion_coef_id, const int grad_q_id, const int face_coef_id) +{ + // Check for negative diffusion coefficients. We don't like them. + // Zero is ok since epsilon^2 is added later + math::HierarchySideDataOpsReal sideops( + hierarchy, 0, -hierarchy->getFinestLevelNumber()); + double diffusion_coef_min = sideops.min(diffusion_coef_id); + if ((diffusion_coef_min + d_epsilon_q * d_epsilon_q) < 0.) { + TBOX_ERROR( + "QuatFaceCoeffs::computeCoeffs(), Negative diffusion coefficient!"); + } + + for (int ln = 0; ln <= hierarchy->getFinestLevelNumber(); ++ln) { + std::shared_ptr level = hierarchy->getPatchLevel(ln); + + for (hier::PatchLevel::Iterator pi(level->begin()); pi != level->end(); + pi++) { + std::shared_ptr patch = *pi; + + std::shared_ptr > diffusion_coef_data( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch->getPatchData(diffusion_coef_id))); + std::shared_ptr > grad_q_data( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch->getPatchData(grad_q_id))); + std::shared_ptr > face_coef_data( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch->getPatchData(face_coef_id))); + + assert(diffusion_coef_data->getDepth() == 1); + computeCoeffs(*patch, *diffusion_coef_data, *grad_q_data, + *face_coef_data); + } + } +} + +// face_coef_data: output +void QuatFaceCoeffs::computeCoeffs(const hier::Patch& patch, + pdat::SideData& diffusion_coef_data, + pdat::SideData& grad_q_data, + pdat::SideData& face_coef_data) +{ + assert(patch.inHierarchy()); + assert(diffusion_coef_data.getDepth() == 1); + assert(grad_q_data.getDepth() == NDIM * d_qlen); + assert(face_coef_data.getDepth() == d_qlen); + + const hier::Box& box = patch.getBox(); + const hier::Index& lower = box.lower(); + const hier::Index& upper = box.upper(); + + const hier::Box& dc_gbox = diffusion_coef_data.getGhostBox(); + const hier::Index& dcglower = dc_gbox.lower(); + const hier::Index& dcgupper = dc_gbox.upper(); + + const hier::Box& gq_gbox = grad_q_data.getGhostBox(); + const hier::Index& gqlower = gq_gbox.lower(); + const hier::Index& gqupper = gq_gbox.upper(); + + const hier::Box& d_gbox = face_coef_data.getGhostBox(); + const hier::Index& dlower = d_gbox.lower(); + const hier::Index& dupper = d_gbox.upper(); + +#if NDIM == 2 + COMPUTE_FACE_COEF2D(lower[0], upper[0], lower[1], upper[1], d_qlen, + d_epsilon_q, diffusion_coef_data.getPointer(0), + dcglower[0], dcgupper[0] + 1, dcglower[1], dcgupper[1], + diffusion_coef_data.getPointer(1), dcglower[0], + dcgupper[0], dcglower[1], dcgupper[1] + 1, + grad_q_data.getPointer(0), gqlower[0], gqupper[0] + 1, + gqlower[1], gqupper[1], grad_q_data.getPointer(1), + gqlower[0], gqupper[0], gqlower[1], gqupper[1] + 1, + face_coef_data.getPointer(0), dlower[0], dupper[0] + 1, + dlower[1], dupper[1], // output + face_coef_data.getPointer(1), dlower[0], dupper[0], + dlower[1], dupper[1] + 1, // output + d_gradient_floor, d_grad_floor_type.c_str()); +#endif +#if NDIM == 3 + COMPUTE_FACE_COEF3D( + lower[0], upper[0], lower[1], upper[1], lower[2], upper[2], d_qlen, + d_epsilon_q, diffusion_coef_data.getPointer(0), dcglower[0], + dcgupper[0] + 1, dcglower[1], dcgupper[1], dcglower[2], dcgupper[2], + diffusion_coef_data.getPointer(1), dcglower[0], dcgupper[0], dcglower[1], + dcgupper[1] + 1, dcglower[2], dcgupper[2], + diffusion_coef_data.getPointer(2), dcglower[0], dcgupper[0], dcglower[1], + dcgupper[1], dcglower[2], dcgupper[2] + 1, grad_q_data.getPointer(0), + gqlower[0], gqupper[0] + 1, gqlower[1], gqupper[1], gqlower[2], + gqupper[2], grad_q_data.getPointer(1), gqlower[0], gqupper[0], + gqlower[1], gqupper[1] + 1, gqlower[2], gqupper[2], + grad_q_data.getPointer(2), gqlower[0], gqupper[0], gqlower[1], + gqupper[1], gqlower[2], gqupper[2] + 1, face_coef_data.getPointer(0), + dlower[0], dupper[0] + 1, dlower[1], dupper[1], dlower[2], dupper[2], + face_coef_data.getPointer(1), dlower[0], dupper[0], dlower[1], + dupper[1] + 1, dlower[2], dupper[2], face_coef_data.getPointer(2), + dlower[0], dupper[0], dlower[1], dupper[1], dlower[2], dupper[2] + 1, + d_gradient_floor, d_grad_floor_type.c_str()); +#endif +} diff --git a/source/QuatFaceCoeffs.h b/source/QuatFaceCoeffs.h new file mode 100644 index 00000000..99a857e0 --- /dev/null +++ b/source/QuatFaceCoeffs.h @@ -0,0 +1,44 @@ +// Copyright (c) 2018, Lawrence Livermore National Security, LLC and +// UT-Battelle, LLC. +// Produced at the Lawrence Livermore National Laboratory and +// the Oak Ridge National Laboratory +// LLNL-CODE-747500 +// All rights reserved. +// This file is part of AMPE. +// For details, see https://github.com/LLNL/AMPE +// Please also read AMPE/LICENSE. +#ifndef included_SAMRAI_config +#include "SAMRAI/SAMRAI_config.h" +#endif + +#include "SAMRAI/hier/Patch.h" +#include "SAMRAI/hier/PatchHierarchy.h" +#include "SAMRAI/pdat/SideData.h" + +#include +#include + +class QuatFaceCoeffs +{ + public: + QuatFaceCoeffs(const int qlen, const double epsilon_q, + const double gradient_floor, + const std::string grad_floor_type); + + // Evaluate coefficient eps_q^2+D_q(phi)/|nabla q| + void computeCoeffs( + const std::shared_ptr hierarchy, + const int diffusion_coef_id, const int grad_q_id, + const int face_coef_id); + + private: + const int d_qlen; + const double d_epsilon_q; + const double d_gradient_floor; + const std::string d_grad_floor_type; + + void computeCoeffs(const SAMRAI::hier::Patch& patch, + SAMRAI::pdat::SideData& diffusion_coef_data, + SAMRAI::pdat::SideData& grad_q_data, + SAMRAI::pdat::SideData& face_coef_data); +}; diff --git a/source/QuatIntegrator.cc b/source/QuatIntegrator.cc index f9ffff58..7a071fd5 100644 --- a/source/QuatIntegrator.cc +++ b/source/QuatIntegrator.cc @@ -142,6 +142,7 @@ QuatIntegrator::QuatIntegrator( d_quat_mobility_id(-1), d_quat_diffusion_id(-1), d_quat_diffusion_deriv_id(-1), + d_quat_face_coef_id(-1), d_quat_symm_rotation_id(-1), d_conc_diffusion_id(-1), d_conc_phase_coupling_diffusion_id(-1), @@ -501,7 +502,9 @@ void QuatIntegrator::setupPreconditioners() d_name + "_QuatIntegratorQuatSy" "sSolver", quatsys_db)); - + d_quat_face_coeffs.reset(new QuatFaceCoeffs(d_qlen, d_epsilon_q, + d_quat_grad_floor, + d_quat_smooth_floor_type)); assert(d_quat_sys_solver); } else { d_quat_sys_solver.reset(); @@ -767,7 +770,6 @@ void QuatIntegrator::RegisterQuatVariables( d_quat_diffs_var, d_current, hier::IntVector(tbox::Dimension(NDIM), NGHOSTS)); - if (d_symmetry_aware) { assert(quat_symm_rotation_var); @@ -1189,8 +1191,17 @@ void QuatIntegrator::RegisterLocalQuatVariables() assert(d_quat_rhs_id >= 0); d_local_data.setFlag(d_quat_rhs_id); - if (d_precond_has_dquatdphi) { + // coeff eps^2+D(\phi)/|grad| to be used in diffusion operator + d_quat_face_coef_var.reset( + new pdat::SideVariable(tbox::Dimension(NDIM), + d_name + "_QI_quat_face_coef_", d_qlen)); + d_quat_face_coef_id = variable_db->registerVariableAndContext( + d_quat_face_coef_var, d_current, + hier::IntVector(tbox::Dimension(NDIM), 0)); + assert(d_quat_face_coef_id >= 0); + d_local_data.setFlag(d_quat_face_coef_id); + if (d_precond_has_dquatdphi) { d_quat_mobility_deriv_var.reset(new pdat::CellVariable( tbox::Dimension(NDIM), d_name + "_QI_quat_mobility_deriv_", 1)); d_quat_mobility_deriv_id = variable_db->registerVariableAndContext( @@ -2783,13 +2794,17 @@ void QuatIntegrator::evaluateQuatRHS( int quat_symm_rotation_id = d_use_gradq_for_flux ? d_quat_symm_rotation_id : -1; + d_quat_face_coeffs->computeCoeffs(hierarchy, d_quat_diffusion_id, + d_quat_grad_side_copy_id, + d_quat_face_coef_id); + // compute RHS using the gradient of q at sides (d_quat_grad_side_id) // computed with physical BC d_quat_sys_solver->evaluateRHS(d_epsilon_q, d_quat_grad_floor, - d_quat_smooth_floor_type, d_quat_diffusion_id, + d_quat_smooth_floor_type, d_quat_face_coef_id, d_quat_grad_side_id, d_quat_grad_side_copy_id, quat_symm_rotation_id, d_quat_mobility_id, - d_quat_scratch_id, quat_rhs_id, true); + d_quat_scratch_id, quat_rhs_id); if (visit_flag && d_model_parameters.with_rhs_visit_output()) @@ -3452,11 +3467,12 @@ int QuatIntegrator::CVSpgmrPrecondSet(double t, SundialsAbstractVector* y, if (d_evolve_quat) { assert(d_quat_sys_solver); - d_quat_sys_solver->setOperatorCoefficients( - gamma, d_epsilon_q, d_quat_grad_floor, d_quat_smooth_floor_type, - d_quat_mobility_id, d_quat_mobility_deriv_id, d_quat_diffusion_id, - d_quat_diffusion_deriv_id, d_quat_grad_side_copy_id, - d_quat_scratch_id); + d_quat_sys_solver->setOperatorCoefficients(gamma, d_quat_mobility_id, + d_quat_mobility_deriv_id, + d_quat_face_coef_id, + d_quat_diffusion_deriv_id, + d_quat_grad_side_copy_id, + d_quat_scratch_id); } // Tell the integrator that the Jacobian data was recomputed @@ -3700,7 +3716,8 @@ int QuatIntegrator::QuatPrecondSolve( // Compute the product of DQuatDPhi block of the Jacobian with the // just computed phi correction - d_quat_sys_solver->multiplyDQuatDPhiBlock(d_phase_sol_id, d_quat_rhs_id); + d_quat_sys_solver->multiplyDQuatDPhiBlock(d_phase_sol_id, d_quat_rhs_id, + d_quat_face_coef_id); // Add gamma times the just computed product to the right-hand side cellops.axpy(d_quat_rhs_id, gamma, d_quat_rhs_id, r_quat_id, false); diff --git a/source/QuatIntegrator.h b/source/QuatIntegrator.h index 0e5b0000..a3ac0a6a 100644 --- a/source/QuatIntegrator.h +++ b/source/QuatIntegrator.h @@ -29,6 +29,7 @@ #include "TemperatureRHSStrategy.h" #include "MovingFrameRHS.h" #include "AdaptMovingFrame.h" +#include "QuatFaceCoeffs.h" // Headers for SAMRAI objects #include "SAMRAI/tbox/Database.h" @@ -677,6 +678,8 @@ class QuatIntegrator : public mesh::StandardTagAndInitStrategy, std::shared_ptr > d_quat_diffs_var; int d_quat_diffs_id; + std::shared_ptr > d_quat_face_coef_var; + int d_quat_face_coef_id; std::shared_ptr > d_f_l_var; int d_f_l_id; @@ -817,6 +820,7 @@ class QuatIntegrator : public mesh::StandardTagAndInitStrategy, CVODESolver* d_sundials_solver; std::shared_ptr d_quat_sys_solver; + std::shared_ptr d_quat_face_coeffs; std::shared_ptr d_phase_fac_ops; std::shared_ptr d_phase_sys_solver; diff --git a/source/QuatLevelSolver.cc b/source/QuatLevelSolver.cc index 40c4b75b..109636f3 100644 --- a/source/QuatLevelSolver.cc +++ b/source/QuatLevelSolver.cc @@ -824,6 +824,7 @@ void QuatLevelSolver::setMatrixCoefficients(const double gamma, std::shared_ptr > face_coef_data( SAMRAI_SHARED_PTR_CAST, hier::PatchData>( patch.getPatchData(face_coef_id))); + assert(face_coef_data->getDepth() == d_qlen); // Set the J_ij blocks for (int depth = 0; depth < d_qlen; depth++) { diff --git a/source/QuatSysSolver.cc b/source/QuatSysSolver.cc index 788f6fe1..49920322 100644 --- a/source/QuatSysSolver.cc +++ b/source/QuatSysSolver.cc @@ -2,7 +2,6 @@ // UT-Battelle, LLC. // Produced at the Lawrence Livermore National Laboratory and // the Oak Ridge National Laboratory -// Written by M.R. Dorr, J.-L. Fattebert and M.E. Wickett // LLNL-CODE-747500 // All rights reserved. // This file is part of AMPE. @@ -264,11 +263,12 @@ void QuatSysSolver::setBoundaries(const std::string& boundary_type, } -void QuatSysSolver::setOperatorCoefficients( - const double time_step, const double epsilon_q, - const double quat_grad_floor, const std::string quat_smooth_floor_type, - const int mobility_id, const int mobility_deriv_id, const int diff_coef_id, - const int diff_coef_deriv_id, const int grad_q_id, const int q_id) +void QuatSysSolver::setOperatorCoefficients(const double time_step, + const int mobility_id, + const int mobility_deriv_id, + const int diff_coef_id, + const int diff_coef_deriv_id, + const int grad_q_id, const int q_id) { t_set_op_coef->start(); @@ -285,10 +285,9 @@ void QuatSysSolver::setOperatorCoefficients( } #endif - d_fac_ops->setOperatorCoefficients(time_step, epsilon_q, mobility_id, - mobility_deriv_id, diff_coef_id, - diff_coef_deriv_id, grad_q_id, q_id, - quat_grad_floor, quat_smooth_floor_type); + d_fac_ops->setOperatorCoefficients(time_step, mobility_id, mobility_deriv_id, + diff_coef_id, diff_coef_deriv_id, + grad_q_id, q_id); t_set_op_coef->stop(); } @@ -356,7 +355,7 @@ void QuatSysSolver::evaluateRHS( const int grad_q_id, const int grad_q_copy_id, // for computation of diffusion coefficient const int rotations_id, const int mobility_id, const int solution_id, - int rhs_id, const bool use_gradq_for_flux) + int rhs_id) { assert(diffusion_coef_id >= 0); assert(grad_q_id >= 0); @@ -366,18 +365,18 @@ void QuatSysSolver::evaluateRHS( d_fac_ops->evaluateRHS(epsilon_q, diffusion_coef_id, grad_q_id, grad_q_copy_id, quat_grad_floor, quat_floor_type, - mobility_id, rotations_id, solution_id, rhs_id, - use_gradq_for_flux); + mobility_id, rotations_id, solution_id, rhs_id); } void QuatSysSolver::multiplyDQuatDPhiBlock(const int q_id, - const int operator_q_id) + const int operator_q_id, + const int face_coef_id) { assert(q_id >= 0); assert(operator_q_id >= 0); - d_fac_ops->multiplyDQuatDPhiBlock(q_id, operator_q_id); + d_fac_ops->multiplyDQuatDPhiBlock(q_id, operator_q_id, face_coef_id); } diff --git a/source/QuatSysSolver.h b/source/QuatSysSolver.h index 194fc834..088769d0 100644 --- a/source/QuatSysSolver.h +++ b/source/QuatSysSolver.h @@ -2,7 +2,6 @@ // UT-Battelle, LLC. // Produced at the Lawrence Livermore National Laboratory and // the Oak Ridge National Laboratory -// Written by M.R. Dorr, J.-L. Fattebert and M.E. Wickett // LLNL-CODE-747500 // All rights reserved. // This file is part of AMPE. @@ -263,10 +262,7 @@ class QuatSysSolver d_fac_ops->setPhysicalBcCoefObject(d_bc_object); } - void setOperatorCoefficients(const double time_step, const double epsilon_q, - const double quat_grad_floor, - const std::string quat_smooth_floor_type, - const int mobility_id, + void setOperatorCoefficients(const double time_step, const int mobility_id, const int mobility_deriv_id, const int diff_coef_id, const int diff_coef_deriv_id, @@ -278,10 +274,10 @@ class QuatSysSolver const std::string quat_smooth_floor_type, const int diff_coef_id, const int grad_q_id, const int grad_q_copy_id, const int rotations_id, - const int mobility_id, const int solution_id, int rhs_id, - const bool use_gradq_for_flux = false); + const int mobility_id, const int solution_id, int rhs_id); - void multiplyDQuatDPhiBlock(const int q_id, const int operator_q_id); + void multiplyDQuatDPhiBlock(const int q_id, const int operator_q_id, + const int face_coef_id); void applyProjection(const int q_id, const int corr_id, const int err_id);