Skip to content

Commit fe0b3fd

Browse files
jeanlucf22jeanlucf22
committed
Compute quaternion face diffusion coeffs outside FACOps
1 parent 04b0cb2 commit fe0b3fd

10 files changed

+240
-180
lines changed

‎source/CMakeLists.txt

+1
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@ set(SOURCES_WPENALTY
2727

2828
set(SOURCES ${CMAKE_SOURCE_DIR}/source/KKStools.cc
2929
${CMAKE_SOURCE_DIR}/source/main.cc
30+
${CMAKE_SOURCE_DIR}/source/QuatFaceCoeffs.cc
3031
${CMAKE_SOURCE_DIR}/source/AdaptMovingFrame.cc
3132
${CMAKE_SOURCE_DIR}/source/computeQDiffs.cc
3233
${CMAKE_SOURCE_DIR}/source/ThreeArgsInterpolationType.cc

‎source/QuatFACOps.cc

+13-148
Original file line numberDiff line numberDiff line change
@@ -680,58 +680,6 @@ void QuatFACOps::initializeOperatorState(
680680
}
681681
}
682682

683-
684-
void QuatFACOps::computeFaceCoefs(const double epsilon_q,
685-
const int diffusion_coef_id,
686-
const int grad_q_id,
687-
const double gradient_floor,
688-
const std::string grad_floor_type,
689-
const int face_coef_id) // output
690-
{
691-
692-
// Check for negative diffusion coefficients. We don't like them.
693-
// Zero is ok since epsilon^2 is added later
694-
695-
math::HierarchySideDataOpsReal<double> sideops(d_hierarchy, d_ln_min,
696-
d_ln_max - 1);
697-
698-
double diffusion_coef_min = sideops.min(diffusion_coef_id);
699-
700-
if ((diffusion_coef_min + epsilon_q * epsilon_q) < 0.) {
701-
TBOX_ERROR(d_object_name << ": Negative diffusion coefficient passed to "
702-
"computeFaceCoefs().");
703-
}
704-
705-
for (int ln = d_ln_min; ln <= d_ln_max; ++ln) {
706-
std::shared_ptr<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);
707-
708-
for (hier::PatchLevel::Iterator pi(level->begin()); pi != level->end();
709-
pi++) {
710-
std::shared_ptr<hier::Patch> patch = *pi;
711-
712-
std::shared_ptr<pdat::SideData<double> > diffusion_coef_data(
713-
SAMRAI_SHARED_PTR_CAST<pdat::SideData<double>, hier::PatchData>(
714-
patch->getPatchData(diffusion_coef_id)));
715-
std::shared_ptr<pdat::SideData<double> > grad_q_data(
716-
SAMRAI_SHARED_PTR_CAST<pdat::SideData<double>, hier::PatchData>(
717-
patch->getPatchData(grad_q_id)));
718-
std::shared_ptr<pdat::SideData<double> > face_coef_data(
719-
SAMRAI_SHARED_PTR_CAST<pdat::SideData<double>, hier::PatchData>(
720-
patch->getPatchData(face_coef_id)));
721-
722-
assert(diffusion_coef_data->getDepth() == 1);
723-
computeFaceCoefsOnPatch(*patch, epsilon_q, *diffusion_coef_data,
724-
*grad_q_data, *face_coef_data, gradient_floor,
725-
grad_floor_type);
726-
}
727-
}
728-
#if 0
729-
double norm = sideops.L1Norm( face_coef_id );
730-
tbox::plog<<"L1 Norm face_coef_id="<<norm<<endl;
731-
#endif
732-
}
733-
734-
735683
void QuatFACOps::computeDQuatDPhiFaceCoefs(const int dprime_id,
736684
const int phi_id,
737685
const int face_coef_id)
@@ -782,14 +730,13 @@ void QuatFACOps::takeSquareRootOnPatch(pdat::CellData<double>& data)
782730
}
783731

784732

785-
void QuatFACOps::setOperatorCoefficients(
786-
const double gamma, const double epsilon_q, const int mobility_id,
787-
const int mobility_deriv_id, const int diffusion_coef_id,
788-
const int face_coef_deriv_id, const int grad_q_id, const int q_id,
789-
const double gradient_floor, const std::string grad_floor_type)
733+
void QuatFACOps::setOperatorCoefficients(const double gamma,
734+
const int mobility_id,
735+
const int mobility_deriv_id,
736+
const int face_coef_id,
737+
const int face_coef_deriv_id,
738+
const int grad_q_id, const int q_id)
790739
{
791-
assert(epsilon_q > 0.);
792-
793740
d_gamma = gamma;
794741

795742
// Check for non-positive mobility
@@ -800,11 +747,6 @@ void QuatFACOps::setOperatorCoefficients(
800747
"setOperatorCoefficients().");
801748
}
802749

803-
// Compute the face coefficients
804-
805-
computeFaceCoefs(epsilon_q, diffusion_coef_id, grad_q_id, gradient_floor,
806-
grad_floor_type, d_face_coef_id);
807-
808750
for (int ln = d_ln_min; ln <= d_ln_max; ++ln) {
809751
std::shared_ptr<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);
810752

@@ -867,7 +809,7 @@ void QuatFACOps::setOperatorCoefficients(
867809
// Set the matrix coefficients
868810
d_quat_level_solver[ln - d_ln_min]->setMatrixCoefficients(d_gamma,
869811
d_sqrt_m_id,
870-
d_face_coef_id);
812+
face_coef_id);
871813
}
872814
}
873815

@@ -1372,81 +1314,7 @@ void QuatFACOps::checkFluxPatchDataIndex() const
13721314
}
13731315
}
13741316

1375-
1376-
/*
1377-
*******************************************************************
1378-
* *
1379-
* AMR-unaware patch-centered computational kernels. *
1380-
* *
1381-
*******************************************************************
1382-
*/
1383-
void QuatFACOps::computeFaceCoefsOnPatch(
1384-
const hier::Patch& patch, const double epsilon_q,
1385-
pdat::SideData<double>& diffusion_coef_data,
1386-
pdat::SideData<double>& grad_q_data,
1387-
pdat::SideData<double>& face_coef_data, // output
1388-
const double gradient_floor, const std::string grad_floor_type) const
1389-
{
1390-
#ifdef DEBUG_CHECK_ASSERTIONS
1391-
assert(patch.inHierarchy());
1392-
assert(diffusion_coef_data.getDepth() == 1);
1393-
assert(grad_q_data.getDepth() == NDIM * d_qlen);
1394-
assert(face_coef_data.getDepth() == d_qlen);
1395-
#endif
1396-
1397-
const hier::Box& box = patch.getBox();
1398-
const hier::Index& lower = box.lower();
1399-
const hier::Index& upper = box.upper();
1400-
1401-
const hier::Box& dc_gbox = diffusion_coef_data.getGhostBox();
1402-
const hier::Index& dcglower = dc_gbox.lower();
1403-
const hier::Index& dcgupper = dc_gbox.upper();
1404-
1405-
const hier::Box& gq_gbox = grad_q_data.getGhostBox();
1406-
const hier::Index& gqlower = gq_gbox.lower();
1407-
const hier::Index& gqupper = gq_gbox.upper();
1408-
1409-
const hier::Box& d_gbox = face_coef_data.getGhostBox();
1410-
const hier::Index& dlower = d_gbox.lower();
1411-
const hier::Index& dupper = d_gbox.upper();
1412-
1413-
#if NDIM == 2
1414-
COMPUTE_FACE_COEF2D(lower[0], upper[0], lower[1], upper[1], d_qlen,
1415-
epsilon_q, diffusion_coef_data.getPointer(0),
1416-
dcglower[0], dcgupper[0] + 1, dcglower[1], dcgupper[1],
1417-
diffusion_coef_data.getPointer(1), dcglower[0],
1418-
dcgupper[0], dcglower[1], dcgupper[1] + 1,
1419-
grad_q_data.getPointer(0), gqlower[0], gqupper[0] + 1,
1420-
gqlower[1], gqupper[1], grad_q_data.getPointer(1),
1421-
gqlower[0], gqupper[0], gqlower[1], gqupper[1] + 1,
1422-
face_coef_data.getPointer(0), dlower[0], dupper[0] + 1,
1423-
dlower[1], dupper[1], // output
1424-
face_coef_data.getPointer(1), dlower[0], dupper[0],
1425-
dlower[1], dupper[1] + 1, // output
1426-
gradient_floor, grad_floor_type.c_str());
1427-
#endif
1428-
#if NDIM == 3
1429-
COMPUTE_FACE_COEF3D(
1430-
lower[0], upper[0], lower[1], upper[1], lower[2], upper[2], d_qlen,
1431-
epsilon_q, diffusion_coef_data.getPointer(0), dcglower[0],
1432-
dcgupper[0] + 1, dcglower[1], dcgupper[1], dcglower[2], dcgupper[2],
1433-
diffusion_coef_data.getPointer(1), dcglower[0], dcgupper[0], dcglower[1],
1434-
dcgupper[1] + 1, dcglower[2], dcgupper[2],
1435-
diffusion_coef_data.getPointer(2), dcglower[0], dcgupper[0], dcglower[1],
1436-
dcgupper[1], dcglower[2], dcgupper[2] + 1, grad_q_data.getPointer(0),
1437-
gqlower[0], gqupper[0] + 1, gqlower[1], gqupper[1], gqlower[2],
1438-
gqupper[2], grad_q_data.getPointer(1), gqlower[0], gqupper[0],
1439-
gqlower[1], gqupper[1] + 1, gqlower[2], gqupper[2],
1440-
grad_q_data.getPointer(2), gqlower[0], gqupper[0], gqlower[1],
1441-
gqupper[1], gqlower[2], gqupper[2] + 1, face_coef_data.getPointer(0),
1442-
dlower[0], dupper[0] + 1, dlower[1], dupper[1], dlower[2], dupper[2],
1443-
face_coef_data.getPointer(1), dlower[0], dupper[0], dlower[1],
1444-
dupper[1] + 1, dlower[2], dupper[2], face_coef_data.getPointer(2),
1445-
dlower[0], dupper[0], dlower[1], dupper[1], dlower[2], dupper[2] + 1,
1446-
gradient_floor, grad_floor_type.c_str());
1447-
#endif
1448-
}
1449-
1317+
// compute face_coef_data using dprime_data and phi_data
14501318
void QuatFACOps::computeDQuatDPhiFaceCoefsOnPatch(
14511319
const hier::Patch& patch, pdat::SideData<double>& dprime_data,
14521320
pdat::CellData<double>& phi_data,
@@ -2009,20 +1877,17 @@ void QuatFACOps::evaluateRHS(
20091877
// Initialize the output array
20101878
d_hopscell->setToScalar(rhs_id, 0., false);
20111879

2012-
computeFaceCoefs(epsilon_q, diffusion_coef_id, grad_q_copy_id,
2013-
gradient_floor, gradient_floor_type,
2014-
d_face_coef_scratch_id);
2015-
20161880
for (int ln = d_ln_max; ln >= d_ln_min; ln--) {
2017-
accumulateOperatorOnLevel(mobility_id, d_face_coef_scratch_id, q_id,
2018-
grad_q_id, rhs_id, ln, true, false);
1881+
accumulateOperatorOnLevel(mobility_id, diffusion_coef_id, q_id, grad_q_id,
1882+
rhs_id, ln, true, false);
20191883
}
20201884

20211885
t_compute_rhs->stop();
20221886
}
20231887

20241888

2025-
void QuatFACOps::multiplyDQuatDPhiBlock(const int phase_id, const int out_id)
1889+
void QuatFACOps::multiplyDQuatDPhiBlock(const int phase_id, const int out_id,
1890+
const int face_coef_id)
20261891
{
20271892
// Initialize the output array
20281893
d_hopscell->setToScalar(out_id, 0., false);
@@ -2031,7 +1896,7 @@ void QuatFACOps::multiplyDQuatDPhiBlock(const int phase_id, const int out_id)
20311896
// mobility derivative times the input phi
20321897
for (int ln = d_ln_max; ln >= d_ln_min; ln--) {
20331898

2034-
accumulateOperatorOnLevel(d_m_deriv_id, d_face_coef_id, d_q_local_id, -1,
1899+
accumulateOperatorOnLevel(d_m_deriv_id, face_coef_id, d_q_local_id, -1,
20351900
out_id, ln, false, false);
20361901

20371902
std::shared_ptr<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);

‎source/QuatFACOps.h

+7-6
Original file line numberDiff line numberDiff line change
@@ -287,11 +287,11 @@ class QuatFACOps : public SAMRAI::solv::FACOperatorStrategy
287287
/*
288288
* Set the operator coefficients.
289289
*/
290-
void setOperatorCoefficients(
291-
const double time_step, const double epsilon_q, const int mobility_id,
292-
const int mobility_deriv_id, const int diff_coef_id,
293-
const int diff_coef_deriv_id, const int grad_q_id, const int q_id,
294-
const double gradient_floor, const std::string grad_floor_type);
290+
void setOperatorCoefficients(const double time_step, const int mobility_id,
291+
const int mobility_deriv_id,
292+
const int diff_coef_id,
293+
const int diff_coef_deriv_id,
294+
const int grad_q_id, const int q_id);
295295

296296
// FACOperatorStrategy virtuals
297297

@@ -353,7 +353,8 @@ class QuatFACOps : public SAMRAI::solv::FACOperatorStrategy
353353
int operator_q_id, int ln, bool project,
354354
bool error_equation_indicator);
355355

356-
void multiplyDQuatDPhiBlock(const int q_id, const int operator_q_id);
356+
void multiplyDQuatDPhiBlock(const int q_id, const int operator_q_id,
357+
const int face_coeff_id);
357358

358359
void applyProjectionOnLevel(const int q_id, const int corr_id,
359360
const int err_id, const int ln);

‎source/QuatFaceCoeffs.cc

+128
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,128 @@
1+
// Copyright (c) 2018, Lawrence Livermore National Security, LLC and
2+
// UT-Battelle, LLC.
3+
// Produced at the Lawrence Livermore National Laboratory and
4+
// the Oak Ridge National Laboratory
5+
// LLNL-CODE-747500
6+
// All rights reserved.
7+
// This file is part of AMPE.
8+
// For details, see https://github.com/LLNL/AMPE
9+
// Please also read AMPE/LICENSE.
10+
#include "QuatFaceCoeffs.h"
11+
#include "QuatFort.h"
12+
13+
#include "SAMRAI/math/HierarchySideDataOpsReal.h"
14+
15+
using namespace SAMRAI;
16+
17+
QuatFaceCoeffs::QuatFaceCoeffs(const int qlen, const double epsilon_q,
18+
const double gradient_floor,
19+
const std::string grad_floor_type)
20+
: d_qlen(qlen),
21+
d_epsilon_q(epsilon_q),
22+
d_gradient_floor(gradient_floor),
23+
d_grad_floor_type(grad_floor_type)
24+
{
25+
}
26+
27+
// Evaluate coefficient eps_q^2+D_q(phi)/|nabla q|
28+
void QuatFaceCoeffs::computeCoeffs(
29+
const std::shared_ptr<hier::PatchHierarchy> hierarchy,
30+
const int diffusion_coef_id, const int grad_q_id, const int face_coef_id)
31+
{
32+
// Check for negative diffusion coefficients. We don't like them.
33+
// Zero is ok since epsilon^2 is added later
34+
math::HierarchySideDataOpsReal<double> sideops(
35+
hierarchy, 0, -hierarchy->getFinestLevelNumber());
36+
double diffusion_coef_min = sideops.min(diffusion_coef_id);
37+
if ((diffusion_coef_min + d_epsilon_q * d_epsilon_q) < 0.) {
38+
TBOX_ERROR(
39+
"QuatFaceCoeffs::computeCoeffs(), Negative diffusion coefficient!");
40+
}
41+
42+
for (int ln = 0; ln <= hierarchy->getFinestLevelNumber(); ++ln) {
43+
std::shared_ptr<hier::PatchLevel> level = hierarchy->getPatchLevel(ln);
44+
45+
for (hier::PatchLevel::Iterator pi(level->begin()); pi != level->end();
46+
pi++) {
47+
std::shared_ptr<hier::Patch> patch = *pi;
48+
49+
std::shared_ptr<pdat::SideData<double> > diffusion_coef_data(
50+
SAMRAI_SHARED_PTR_CAST<pdat::SideData<double>, hier::PatchData>(
51+
patch->getPatchData(diffusion_coef_id)));
52+
std::shared_ptr<pdat::SideData<double> > grad_q_data(
53+
SAMRAI_SHARED_PTR_CAST<pdat::SideData<double>, hier::PatchData>(
54+
patch->getPatchData(grad_q_id)));
55+
std::shared_ptr<pdat::SideData<double> > face_coef_data(
56+
SAMRAI_SHARED_PTR_CAST<pdat::SideData<double>, hier::PatchData>(
57+
patch->getPatchData(face_coef_id)));
58+
59+
assert(diffusion_coef_data->getDepth() == 1);
60+
computeCoeffs(*patch, *diffusion_coef_data, *grad_q_data,
61+
*face_coef_data);
62+
}
63+
}
64+
}
65+
66+
// face_coef_data: output
67+
void QuatFaceCoeffs::computeCoeffs(const hier::Patch& patch,
68+
pdat::SideData<double>& diffusion_coef_data,
69+
pdat::SideData<double>& grad_q_data,
70+
pdat::SideData<double>& face_coef_data)
71+
{
72+
assert(patch.inHierarchy());
73+
assert(diffusion_coef_data.getDepth() == 1);
74+
assert(grad_q_data.getDepth() == NDIM * d_qlen);
75+
assert(face_coef_data.getDepth() == d_qlen);
76+
77+
const hier::Box& box = patch.getBox();
78+
const hier::Index& lower = box.lower();
79+
const hier::Index& upper = box.upper();
80+
81+
const hier::Box& dc_gbox = diffusion_coef_data.getGhostBox();
82+
const hier::Index& dcglower = dc_gbox.lower();
83+
const hier::Index& dcgupper = dc_gbox.upper();
84+
85+
const hier::Box& gq_gbox = grad_q_data.getGhostBox();
86+
const hier::Index& gqlower = gq_gbox.lower();
87+
const hier::Index& gqupper = gq_gbox.upper();
88+
89+
const hier::Box& d_gbox = face_coef_data.getGhostBox();
90+
const hier::Index& dlower = d_gbox.lower();
91+
const hier::Index& dupper = d_gbox.upper();
92+
93+
#if NDIM == 2
94+
COMPUTE_FACE_COEF2D(lower[0], upper[0], lower[1], upper[1], d_qlen,
95+
d_epsilon_q, diffusion_coef_data.getPointer(0),
96+
dcglower[0], dcgupper[0] + 1, dcglower[1], dcgupper[1],
97+
diffusion_coef_data.getPointer(1), dcglower[0],
98+
dcgupper[0], dcglower[1], dcgupper[1] + 1,
99+
grad_q_data.getPointer(0), gqlower[0], gqupper[0] + 1,
100+
gqlower[1], gqupper[1], grad_q_data.getPointer(1),
101+
gqlower[0], gqupper[0], gqlower[1], gqupper[1] + 1,
102+
face_coef_data.getPointer(0), dlower[0], dupper[0] + 1,
103+
dlower[1], dupper[1], // output
104+
face_coef_data.getPointer(1), dlower[0], dupper[0],
105+
dlower[1], dupper[1] + 1, // output
106+
d_gradient_floor, d_grad_floor_type.c_str());
107+
#endif
108+
#if NDIM == 3
109+
COMPUTE_FACE_COEF3D(
110+
lower[0], upper[0], lower[1], upper[1], lower[2], upper[2], d_qlen,
111+
d_epsilon_q, diffusion_coef_data.getPointer(0), dcglower[0],
112+
dcgupper[0] + 1, dcglower[1], dcgupper[1], dcglower[2], dcgupper[2],
113+
diffusion_coef_data.getPointer(1), dcglower[0], dcgupper[0], dcglower[1],
114+
dcgupper[1] + 1, dcglower[2], dcgupper[2],
115+
diffusion_coef_data.getPointer(2), dcglower[0], dcgupper[0], dcglower[1],
116+
dcgupper[1], dcglower[2], dcgupper[2] + 1, grad_q_data.getPointer(0),
117+
gqlower[0], gqupper[0] + 1, gqlower[1], gqupper[1], gqlower[2],
118+
gqupper[2], grad_q_data.getPointer(1), gqlower[0], gqupper[0],
119+
gqlower[1], gqupper[1] + 1, gqlower[2], gqupper[2],
120+
grad_q_data.getPointer(2), gqlower[0], gqupper[0], gqlower[1],
121+
gqupper[1], gqlower[2], gqupper[2] + 1, face_coef_data.getPointer(0),
122+
dlower[0], dupper[0] + 1, dlower[1], dupper[1], dlower[2], dupper[2],
123+
face_coef_data.getPointer(1), dlower[0], dupper[0], dlower[1],
124+
dupper[1] + 1, dlower[2], dupper[2], face_coef_data.getPointer(2),
125+
dlower[0], dupper[0], dlower[1], dupper[1], dlower[2], dupper[2] + 1,
126+
d_gradient_floor, d_grad_floor_type.c_str());
127+
#endif
128+
}

0 commit comments

Comments
 (0)