Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Always pass in grad q when computing RHS for q #128

Draft
wants to merge 2 commits into
base: release
Choose a base branch
from
Draft

Always pass in grad q when computing RHS for q #128

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
04b0cb2
Always pass in grad q when computing RHS for q
jeanlucf22 Mar 22, 2023
fe0b3fd
Compute quaternion face diffusion coeffs outside FACOps
jeanlucf22 Mar 22, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions source/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
168 changes: 15 additions & 153 deletions source/QuatFACOps.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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<double> 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<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);

for (hier::PatchLevel::Iterator pi(level->begin()); pi != level->end();
pi++) {
std::shared_ptr<hier::Patch> patch = *pi;

std::shared_ptr<pdat::SideData<double> > diffusion_coef_data(
SAMRAI_SHARED_PTR_CAST<pdat::SideData<double>, hier::PatchData>(
patch->getPatchData(diffusion_coef_id)));
std::shared_ptr<pdat::SideData<double> > grad_q_data(
SAMRAI_SHARED_PTR_CAST<pdat::SideData<double>, hier::PatchData>(
patch->getPatchData(grad_q_id)));
std::shared_ptr<pdat::SideData<double> > face_coef_data(
SAMRAI_SHARED_PTR_CAST<pdat::SideData<double>, 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="<<norm<<endl;
#endif
}


void QuatFACOps::computeDQuatDPhiFaceCoefs(const int dprime_id,
const int phi_id,
const int face_coef_id)
Expand Down Expand Up @@ -783,14 +730,13 @@ void QuatFACOps::takeSquareRootOnPatch(pdat::CellData<double>& 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
Expand All @@ -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<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);

Expand Down Expand Up @@ -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);
}
}

Expand Down Expand Up @@ -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<double>& diffusion_coef_data,
pdat::SideData<double>& grad_q_data,
pdat::SideData<double>& 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<double>& dprime_data,
pdat::CellData<double>& phi_data,
Expand Down Expand Up @@ -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);
Expand All @@ -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<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);
Expand Down
40 changes: 8 additions & 32 deletions source/QuatFACOps.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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);
Expand Down
Loading