Skip to content

Commit 46190ef

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

9 files changed

+202
-18
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

+6-9
Original file line numberDiff line numberDiff line change
@@ -1446,7 +1446,7 @@ void QuatFACOps::computeFaceCoefsOnPatch(
14461446
gradient_floor, grad_floor_type.c_str());
14471447
#endif
14481448
}
1449-
1449+
// compute face_coef_data using dprime_data and phi_data
14501450
void QuatFACOps::computeDQuatDPhiFaceCoefsOnPatch(
14511451
const hier::Patch& patch, pdat::SideData<double>& dprime_data,
14521452
pdat::CellData<double>& phi_data,
@@ -2009,20 +2009,17 @@ void QuatFACOps::evaluateRHS(
20092009
// Initialize the output array
20102010
d_hopscell->setToScalar(rhs_id, 0., false);
20112011

2012-
computeFaceCoefs(epsilon_q, diffusion_coef_id, grad_q_copy_id,
2013-
gradient_floor, gradient_floor_type,
2014-
d_face_coef_scratch_id);
2015-
20162012
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);
2013+
accumulateOperatorOnLevel(mobility_id, diffusion_coef_id, q_id, grad_q_id,
2014+
rhs_id, ln, true, false);
20192015
}
20202016

20212017
t_compute_rhs->stop();
20222018
}
20232019

20242020

2025-
void QuatFACOps::multiplyDQuatDPhiBlock(const int phase_id, const int out_id)
2021+
void QuatFACOps::multiplyDQuatDPhiBlock(const int phase_id, const int out_id,
2022+
const int face_coef_id)
20262023
{
20272024
// Initialize the output array
20282025
d_hopscell->setToScalar(out_id, 0., false);
@@ -2031,7 +2028,7 @@ void QuatFACOps::multiplyDQuatDPhiBlock(const int phase_id, const int out_id)
20312028
// mobility derivative times the input phi
20322029
for (int ln = d_ln_max; ln >= d_ln_min; ln--) {
20332030

2034-
accumulateOperatorOnLevel(d_m_deriv_id, d_face_coef_id, d_q_local_id, -1,
2031+
accumulateOperatorOnLevel(d_m_deriv_id, face_coef_id, d_q_local_id, -1,
20352032
out_id, ln, false, false);
20362033

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

‎source/QuatFACOps.h

+2-1
Original file line numberDiff line numberDiff line change
@@ -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

+119
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,119 @@
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+
14+
using namespace SAMRAI;
15+
16+
QuatFaceCoeffs::QuatFaceCoeffs(const int qlen, const double epsilon_q,
17+
const double gradient_floor,
18+
const std::string grad_floor_type)
19+
: d_qlen(qlen),
20+
d_epsilon_q(epsilon_q),
21+
d_gradient_floor(gradient_floor),
22+
d_grad_floor_type(grad_floor_type)
23+
{
24+
}
25+
26+
// Evaluate coefficient eps_q^2+D_q(phi)/|nabla q|
27+
void QuatFaceCoeffs::computeCoeffs(
28+
const std::shared_ptr<hier::PatchHierarchy> hierarchy,
29+
const int diffusion_coef_id, const int grad_q_id, const int face_coef_id)
30+
{
31+
for (int ln = 0; ln <= hierarchy->getFinestLevelNumber(); ++ln) {
32+
std::shared_ptr<hier::PatchLevel> level = hierarchy->getPatchLevel(ln);
33+
34+
for (hier::PatchLevel::Iterator pi(level->begin()); pi != level->end();
35+
pi++) {
36+
std::shared_ptr<hier::Patch> patch = *pi;
37+
38+
std::shared_ptr<pdat::SideData<double> > diffusion_coef_data(
39+
SAMRAI_SHARED_PTR_CAST<pdat::SideData<double>, hier::PatchData>(
40+
patch->getPatchData(diffusion_coef_id)));
41+
std::shared_ptr<pdat::SideData<double> > grad_q_data(
42+
SAMRAI_SHARED_PTR_CAST<pdat::SideData<double>, hier::PatchData>(
43+
patch->getPatchData(grad_q_id)));
44+
std::shared_ptr<pdat::SideData<double> > face_coef_data(
45+
SAMRAI_SHARED_PTR_CAST<pdat::SideData<double>, hier::PatchData>(
46+
patch->getPatchData(face_coef_id)));
47+
48+
assert(diffusion_coef_data->getDepth() == 1);
49+
computeCoeffs(*patch, *diffusion_coef_data, *grad_q_data,
50+
*face_coef_data);
51+
}
52+
}
53+
}
54+
55+
// face_coef_data: output
56+
void QuatFaceCoeffs::computeCoeffs(const hier::Patch& patch,
57+
pdat::SideData<double>& diffusion_coef_data,
58+
pdat::SideData<double>& grad_q_data,
59+
pdat::SideData<double>& face_coef_data)
60+
{
61+
#ifdef DEBUG_CHECK_ASSERTIONS
62+
assert(patch.inHierarchy());
63+
assert(diffusion_coef_data.getDepth() == 1);
64+
assert(grad_q_data.getDepth() == NDIM * d_qlen);
65+
assert(face_coef_data.getDepth() == d_qlen);
66+
#endif
67+
68+
const hier::Box& box = patch.getBox();
69+
const hier::Index& lower = box.lower();
70+
const hier::Index& upper = box.upper();
71+
72+
const hier::Box& dc_gbox = diffusion_coef_data.getGhostBox();
73+
const hier::Index& dcglower = dc_gbox.lower();
74+
const hier::Index& dcgupper = dc_gbox.upper();
75+
76+
const hier::Box& gq_gbox = grad_q_data.getGhostBox();
77+
const hier::Index& gqlower = gq_gbox.lower();
78+
const hier::Index& gqupper = gq_gbox.upper();
79+
80+
const hier::Box& d_gbox = face_coef_data.getGhostBox();
81+
const hier::Index& dlower = d_gbox.lower();
82+
const hier::Index& dupper = d_gbox.upper();
83+
84+
#if NDIM == 2
85+
COMPUTE_FACE_COEF2D(lower[0], upper[0], lower[1], upper[1], d_qlen,
86+
d_epsilon_q, diffusion_coef_data.getPointer(0),
87+
dcglower[0], dcgupper[0] + 1, dcglower[1], dcgupper[1],
88+
diffusion_coef_data.getPointer(1), dcglower[0],
89+
dcgupper[0], dcglower[1], dcgupper[1] + 1,
90+
grad_q_data.getPointer(0), gqlower[0], gqupper[0] + 1,
91+
gqlower[1], gqupper[1], grad_q_data.getPointer(1),
92+
gqlower[0], gqupper[0], gqlower[1], gqupper[1] + 1,
93+
face_coef_data.getPointer(0), dlower[0], dupper[0] + 1,
94+
dlower[1], dupper[1], // output
95+
face_coef_data.getPointer(1), dlower[0], dupper[0],
96+
dlower[1], dupper[1] + 1, // output
97+
d_gradient_floor, d_grad_floor_type.c_str());
98+
#endif
99+
#if NDIM == 3
100+
COMPUTE_FACE_COEF3D(
101+
lower[0], upper[0], lower[1], upper[1], lower[2], upper[2], d_qlen,
102+
d_epsilon_q, diffusion_coef_data.getPointer(0), dcglower[0],
103+
dcgupper[0] + 1, dcglower[1], dcgupper[1], dcglower[2], dcgupper[2],
104+
diffusion_coef_data.getPointer(1), dcglower[0], dcgupper[0], dcglower[1],
105+
dcgupper[1] + 1, dcglower[2], dcgupper[2],
106+
diffusion_coef_data.getPointer(2), dcglower[0], dcgupper[0], dcglower[1],
107+
dcgupper[1], dcglower[2], dcgupper[2] + 1, grad_q_data.getPointer(0),
108+
gqlower[0], gqupper[0] + 1, gqlower[1], gqupper[1], gqlower[2],
109+
gqupper[2], grad_q_data.getPointer(1), gqlower[0], gqupper[0],
110+
gqlower[1], gqupper[1] + 1, gqlower[2], gqupper[2],
111+
grad_q_data.getPointer(2), gqlower[0], gqupper[0], gqlower[1],
112+
gqupper[1], gqlower[2], gqupper[2] + 1, face_coef_data.getPointer(0),
113+
dlower[0], dupper[0] + 1, dlower[1], dupper[1], dlower[2], dupper[2],
114+
face_coef_data.getPointer(1), dlower[0], dupper[0], dlower[1],
115+
dupper[1] + 1, dlower[2], dupper[2], face_coef_data.getPointer(2),
116+
dlower[0], dupper[0], dlower[1], dupper[1], dlower[2], dupper[2] + 1,
117+
d_gradient_floor, d_grad_floor_type.c_str());
118+
#endif
119+
}

‎source/QuatFaceCoeffs.h

+44
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
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+
#ifndef included_SAMRAI_config
11+
#include "SAMRAI/SAMRAI_config.h"
12+
#endif
13+
14+
#include "SAMRAI/hier/Patch.h"
15+
#include "SAMRAI/hier/PatchHierarchy.h"
16+
#include "SAMRAI/pdat/SideData.h"
17+
18+
#include <string>
19+
#include <memory>
20+
21+
class QuatFaceCoeffs
22+
{
23+
public:
24+
QuatFaceCoeffs(const int qlen, const double epsilon_q,
25+
const double gradient_floor,
26+
const std::string grad_floor_type);
27+
28+
// Evaluate coefficient eps_q^2+D_q(phi)/|nabla q|
29+
void computeCoeffs(
30+
const std::shared_ptr<SAMRAI::hier::PatchHierarchy> hierarchy,
31+
const int diffusion_coef_id, const int grad_q_id,
32+
const int face_coef_id);
33+
34+
private:
35+
const int d_qlen;
36+
const double d_epsilon_q;
37+
const double d_gradient_floor;
38+
const std::string d_grad_floor_type;
39+
40+
void computeCoeffs(const SAMRAI::hier::Patch& patch,
41+
SAMRAI::pdat::SideData<double>& diffusion_coef_data,
42+
SAMRAI::pdat::SideData<double>& grad_q_data,
43+
SAMRAI::pdat::SideData<double>& face_coef_data);
44+
};

‎source/QuatIntegrator.cc

+21-5
Original file line numberDiff line numberDiff line change
@@ -142,6 +142,7 @@ QuatIntegrator::QuatIntegrator(
142142
d_quat_mobility_id(-1),
143143
d_quat_diffusion_id(-1),
144144
d_quat_diffusion_deriv_id(-1),
145+
d_quat_face_coef_id(-1),
145146
d_quat_symm_rotation_id(-1),
146147
d_conc_diffusion_id(-1),
147148
d_conc_phase_coupling_diffusion_id(-1),
@@ -501,7 +502,9 @@ void QuatIntegrator::setupPreconditioners()
501502
d_name + "_QuatIntegratorQuatSy"
502503
"sSolver",
503504
quatsys_db));
504-
505+
d_quat_face_coeffs.reset(new QuatFaceCoeffs(d_qlen, d_epsilon_q,
506+
d_quat_grad_floor,
507+
d_quat_smooth_floor_type));
505508
assert(d_quat_sys_solver);
506509
} else {
507510
d_quat_sys_solver.reset();
@@ -767,7 +770,6 @@ void QuatIntegrator::RegisterQuatVariables(
767770
d_quat_diffs_var, d_current,
768771
hier::IntVector(tbox::Dimension(NDIM), NGHOSTS));
769772

770-
771773
if (d_symmetry_aware) {
772774
assert(quat_symm_rotation_var);
773775

@@ -1189,8 +1191,17 @@ void QuatIntegrator::RegisterLocalQuatVariables()
11891191
assert(d_quat_rhs_id >= 0);
11901192
d_local_data.setFlag(d_quat_rhs_id);
11911193

1192-
if (d_precond_has_dquatdphi) {
1194+
// coeff eps^2+D(\phi)/|grad| to be used in diffusion operator
1195+
d_quat_face_coef_var.reset(
1196+
new pdat::SideVariable<double>(tbox::Dimension(NDIM),
1197+
d_name + "_QI_quat_face_coef_", d_qlen));
1198+
d_quat_face_coef_id = variable_db->registerVariableAndContext(
1199+
d_quat_face_coef_var, d_current,
1200+
hier::IntVector(tbox::Dimension(NDIM), 0));
1201+
assert(d_quat_face_coef_id >= 0);
1202+
d_local_data.setFlag(d_quat_face_coef_id);
11931203

1204+
if (d_precond_has_dquatdphi) {
11941205
d_quat_mobility_deriv_var.reset(new pdat::CellVariable<double>(
11951206
tbox::Dimension(NDIM), d_name + "_QI_quat_mobility_deriv_", 1));
11961207
d_quat_mobility_deriv_id = variable_db->registerVariableAndContext(
@@ -2783,10 +2794,14 @@ void QuatIntegrator::evaluateQuatRHS(
27832794
int quat_symm_rotation_id =
27842795
d_use_gradq_for_flux ? d_quat_symm_rotation_id : -1;
27852796

2797+
d_quat_face_coeffs->computeCoeffs(hierarchy, d_quat_diffusion_id,
2798+
d_quat_grad_side_copy_id,
2799+
d_quat_face_coef_id);
2800+
27862801
// compute RHS using the gradient of q at sides (d_quat_grad_side_id)
27872802
// computed with physical BC
27882803
d_quat_sys_solver->evaluateRHS(d_epsilon_q, d_quat_grad_floor,
2789-
d_quat_smooth_floor_type, d_quat_diffusion_id,
2804+
d_quat_smooth_floor_type, d_quat_face_coef_id,
27902805
d_quat_grad_side_id, d_quat_grad_side_copy_id,
27912806
quat_symm_rotation_id, d_quat_mobility_id,
27922807
d_quat_scratch_id, quat_rhs_id);
@@ -3700,7 +3715,8 @@ int QuatIntegrator::QuatPrecondSolve(
37003715

37013716
// Compute the product of DQuatDPhi block of the Jacobian with the
37023717
// just computed phi correction
3703-
d_quat_sys_solver->multiplyDQuatDPhiBlock(d_phase_sol_id, d_quat_rhs_id);
3718+
d_quat_sys_solver->multiplyDQuatDPhiBlock(d_phase_sol_id, d_quat_rhs_id,
3719+
d_quat_face_coef_id);
37043720

37053721
// Add gamma times the just computed product to the right-hand side
37063722
cellops.axpy(d_quat_rhs_id, gamma, d_quat_rhs_id, r_quat_id, false);

‎source/QuatIntegrator.h

+4
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@
2929
#include "TemperatureRHSStrategy.h"
3030
#include "MovingFrameRHS.h"
3131
#include "AdaptMovingFrame.h"
32+
#include "QuatFaceCoeffs.h"
3233

3334
// Headers for SAMRAI objects
3435
#include "SAMRAI/tbox/Database.h"
@@ -677,6 +678,8 @@ class QuatIntegrator : public mesh::StandardTagAndInitStrategy,
677678

678679
std::shared_ptr<pdat::SideVariable<double> > d_quat_diffs_var;
679680
int d_quat_diffs_id;
681+
std::shared_ptr<pdat::SideVariable<double> > d_quat_face_coef_var;
682+
int d_quat_face_coef_id;
680683

681684
std::shared_ptr<pdat::CellVariable<double> > d_f_l_var;
682685
int d_f_l_id;
@@ -817,6 +820,7 @@ class QuatIntegrator : public mesh::StandardTagAndInitStrategy,
817820
CVODESolver* d_sundials_solver;
818821

819822
std::shared_ptr<QuatSysSolver> d_quat_sys_solver;
823+
std::shared_ptr<QuatFaceCoeffs> d_quat_face_coeffs;
820824

821825
std::shared_ptr<PhaseFACOps> d_phase_fac_ops;
822826
std::shared_ptr<PhaseFACSolver> d_phase_sys_solver;

‎source/QuatSysSolver.cc

+3-2
Original file line numberDiff line numberDiff line change
@@ -370,12 +370,13 @@ void QuatSysSolver::evaluateRHS(
370370

371371

372372
void QuatSysSolver::multiplyDQuatDPhiBlock(const int q_id,
373-
const int operator_q_id)
373+
const int operator_q_id,
374+
const int face_coef_id)
374375
{
375376
assert(q_id >= 0);
376377
assert(operator_q_id >= 0);
377378

378-
d_fac_ops->multiplyDQuatDPhiBlock(q_id, operator_q_id);
379+
d_fac_ops->multiplyDQuatDPhiBlock(q_id, operator_q_id, face_coef_id);
379380
}
380381

381382

‎source/QuatSysSolver.h

+2-1
Original file line numberDiff line numberDiff line change
@@ -279,7 +279,8 @@ class QuatSysSolver
279279
const int grad_q_copy_id, const int rotations_id,
280280
const int mobility_id, const int solution_id, int rhs_id);
281281

282-
void multiplyDQuatDPhiBlock(const int q_id, const int operator_q_id);
282+
void multiplyDQuatDPhiBlock(const int q_id, const int operator_q_id,
283+
const int face_coef_id);
283284

284285
void applyProjection(const int q_id, const int corr_id, const int err_id);
285286

0 commit comments

Comments
 (0)