18#define DefaultBCDirichlet BoundaryCondition{BCType::DIRICHLET,{0,0,0}}
24 Assemble_b(
const std::vector<double>& q_vector)
26 const std::string fname =
"lbs::acceleration::DiffusionMIPSolver::"
28 if (
A_ ==
nullptr or
rhs_ ==
nullptr or
ksp_ ==
nullptr)
29 throw std::logic_error(fname +
": Some or all PETSc elements are null. "
30 "Check that Initialize has been called.");
39 const size_t num_faces = cell.faces_.size();
41 const size_t num_nodes = cell_mapping.
NumNodes();
42 const auto cc_nodes = cell_mapping.GetNodeLocations();
45 const auto& cell_M_matrix = unit_cell_matrices.M_matrix;
49 for (
size_t g=0; g<num_groups; ++g)
52 const double Dg = xs.Dg[g];
54 std::vector<double> qg(num_nodes, 0.0);
55 for (
size_t j=0; j<num_nodes; j++)
59 for (
size_t i=0; i<num_nodes; i++)
62 double entry_rhs_i = 0.0;
63 for (
size_t j=0; j<num_nodes; j++)
64 entry_rhs_i += qg[j]* cell_M_matrix[i][j];
66 VecSetValue(
rhs_, imap, entry_rhs_i, ADD_VALUES);
70 for (
size_t f=0; f<num_faces; ++f)
72 const auto& face = cell.faces_[f];
73 const auto& n_f = face.normal_;
74 const size_t num_face_nodes = cell_mapping.NumFaceNodes(f);
76 const auto& face_M = unit_cell_matrices.face_M_matrices[f];
77 const auto& face_G = unit_cell_matrices.face_G_matrices[f];
78 const auto& face_Si = unit_cell_matrices.face_Si_vectors[f];
82 if (not face.has_neighbor_)
85 if (
bcs_.count(face.neighbor_id_) > 0)
86 bc =
bcs_.at(face.neighbor_id_);
90 const double bc_value = bc.values[0];
102 for (
size_t fi=0; fi<num_face_nodes; ++fi)
104 const int i = cell_mapping.MapFaceNode(f,fi);
107 for (
size_t fj=0; fj<num_face_nodes; ++fj)
109 const int jm = cell_mapping.MapFaceNode(f,fj);
111 const double aij = kappa*face_M[i][jm];
112 const double aij_bc_value = aij*bc_value;
114 VecSetValue(
rhs_, imap, aij_bc_value, ADD_VALUES);
123 for (
size_t i=0; i<num_nodes; i++)
127 for (
size_t j=0; j<num_nodes; j++)
129 const double aij = -Dg*n_f.Dot(face_G[j][i] + face_G[i][j]);
130 const double aij_bc_value = aij*bc_value;
132 VecSetValue(
rhs_, imap, aij_bc_value, ADD_VALUES);
138 const double bval = bc.values[1];
139 const double fval = bc.values[2];
141 if (std::fabs(bval) < 1.0e-12)
continue;
143 for (
size_t fi=0; fi<num_face_nodes; fi++)
145 const int i = cell_mapping.MapFaceNode(f,fi);
148 if (std::fabs(fval) >= 1.0e-12)
150 const double rhs_val = (fval/bval) * face_Si[i];
152 VecSetValue(
rhs_, ir, rhs_val, ADD_VALUES);
161 VecAssemblyBegin(
rhs_);
162 VecAssemblyEnd(
rhs_);
174 const std::string fname =
"lbs::acceleration::DiffusionMIPSolver::"
176 if (A_ ==
nullptr or rhs_ ==
nullptr or ksp_ ==
nullptr)
177 throw std::logic_error(fname +
": Some or all PETSc elements are null. "
178 "Check that Initialize has been called.");
182 const size_t num_groups = uk_man_.unknowns_.front().num_components_;
184 const double* q_vector;
185 VecGetArrayRead(petsc_q_vector, &q_vector);
188 for (
const auto& cell : grid_.local_cells)
190 const size_t num_faces = cell.faces_.size();
191 const auto& cell_mapping = sdm_.GetCellMapping(cell);
192 const size_t num_nodes = cell_mapping.NumNodes();
193 const auto cc_nodes = cell_mapping.GetNodeLocations();
194 const auto& unit_cell_matrices = unit_cell_matrices_[cell.local_id_];
196 const auto& cell_M_matrix = unit_cell_matrices.M_matrix;
198 const auto& xs = mat_id_2_xs_map_.at(cell.material_id_);
200 for (
size_t g=0; g<num_groups; ++g)
203 const double Dg = xs.Dg[g];
205 std::vector<double> qg(num_nodes, 0.0);
206 for (
size_t j=0; j<num_nodes; j++)
207 qg[j] = q_vector[sdm_.MapDOFLocal(cell, j, uk_man_, 0, g)];
210 for (
size_t i=0; i<num_nodes; i++)
212 const int64_t imap = sdm_.MapDOF(cell, i, uk_man_, 0, g);
213 double entry_rhs_i = 0.0;
214 for (
size_t j=0; j<num_nodes; j++)
215 entry_rhs_i += qg[j]* cell_M_matrix[i][j];
217 VecSetValue(rhs_, imap, entry_rhs_i, ADD_VALUES);
221 for (
size_t f=0; f<num_faces; ++f)
223 const auto& face = cell.faces_[f];
224 const auto& n_f = face.normal_;
225 const size_t num_face_nodes = cell_mapping.NumFaceNodes(f);
227 const auto& face_M = unit_cell_matrices.face_M_matrices[f];
228 const auto& face_G = unit_cell_matrices.face_G_matrices[f];
229 const auto& face_Si = unit_cell_matrices.face_Si_vectors[f];
231 const double hm = HPerpendicular(cell, f);
233 if (not face.has_neighbor_)
236 if (bcs_.count(face.neighbor_id_) > 0)
237 bc = bcs_.at(face.neighbor_id_);
241 const double bc_value = bc.values[0];
246 kappa = fmax(options.penalty_factor*Dg/hm,0.25);
248 kappa = fmax(options.penalty_factor*Dg/hm,0.25);
250 kappa = fmax(options.penalty_factor*2.0*Dg/hm,0.25);
253 for (
size_t fi=0; fi<num_face_nodes; ++fi)
255 const int i = cell_mapping.MapFaceNode(f,fi);
256 const int64_t imap = sdm_.MapDOF(cell, i, uk_man_, 0, g);
258 for (
size_t fj=0; fj<num_face_nodes; ++fj)
260 const int jm = cell_mapping.MapFaceNode(f,fj);
262 const double aij = kappa*face_M[i][jm];
263 const double aij_bc_value = aij*bc_value;
265 VecSetValue(rhs_, imap, aij_bc_value, ADD_VALUES);
274 for (
size_t i=0; i<num_nodes; i++)
276 const int64_t imap = sdm_.MapDOF(cell, i, uk_man_, 0, g);
278 for (
size_t j=0; j<num_nodes; j++)
280 const double aij = -Dg*n_f.Dot(face_G[j][i] + face_G[i][j]);
281 const double aij_bc_value = aij*bc_value;
283 VecSetValue(rhs_, imap, aij_bc_value, ADD_VALUES);
289 const double bval = bc.values[1];
290 const double fval = bc.values[2];
292 if (std::fabs(bval) < 1.0e-12)
continue;
294 for (
size_t fi=0; fi<num_face_nodes; fi++)
296 const int i = cell_mapping.MapFaceNode(f,fi);
297 const int64_t ir = sdm_.MapDOF(cell, i, uk_man_, 0, g);
299 if (std::fabs(fval) >= 1.0e-12)
301 const double rhs_val = (fval/bval) * face_Si[i];
303 VecSetValue(rhs_, ir, rhs_val, ADD_VALUES);
312 VecRestoreArrayRead(petsc_q_vector, &q_vector);
314 VecAssemblyBegin(rhs_);
315 VecAssemblyEnd(rhs_);
static chi::Timer program_timer
LogStream Log(LOG_LVL level=LOG_0)
std::string GetTimeString() const
const CellMapping & GetCellMapping(const chi_mesh::Cell &cell) const
virtual int64_t MapDOF(const chi_mesh::Cell &cell, unsigned int node, const UnknownManager &unknown_manager, unsigned int unknown_id, unsigned int component) const =0
virtual int64_t MapDOFLocal(const chi_mesh::Cell &cell, unsigned int node, const UnknownManager &unknown_manager, unsigned int unknown_id, unsigned int component) const =0
std::vector< Unknown > unknowns_
LocalCellHandler local_cells
void Assemble_b(const std::vector< double > &q_vector) override
double HPerpendicular(const chi_mesh::Cell &cell, unsigned int f)
struct lbs::acceleration::DiffusionSolver::Options options
const chi_math::SpatialDiscretization & sdm_
const chi_math::UnknownManager uk_man_
const chi_mesh::MeshContinuum & grid_
const std::map< uint64_t, BoundaryCondition > bcs_
const MatID2XSMap mat_id_2_xs_map_
const std::vector< UnitCellMatrices > & unit_cell_matrices_
#define DefaultBCDirichlet
bool verbose
Verbosity flag.