phynexis-deps/mfem
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'master' into block-solvers/bramble-pasciak-dev

Browse Source
This commit is contained in:
Tzanio Kolev
2024-02-03 14:31:32 -08:00
committed by GitHub
parent 91b9d9dc92 33617ce3a7
commit d53f982761
22 changed files with 388 additions and 75 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fem/bilinearform.cpp
+3 -2
View File
@@ -1585,9 +1585,10 @@ void MixedBilinearForm::Assemble(int skip_zeros)
ftr = mesh->GetBdrFaceTransformations(i);
if (ftr)
{
trial_fes->GetFaceVDofs(ftr->ElementNo, trial_vdofs);
const int iface = mesh->GetBdrElementFaceIndex(i);
trial_fes->GetFaceVDofs(iface, trial_vdofs);
test_fes->GetElementVDofs(ftr->Elem1No, test_vdofs);
trial_face_fe = trial_fes->GetFaceElement(ftr->ElementNo);
trial_face_fe = trial_fes->GetFaceElement(iface);
test_fe1 = test_fes->GetFE(ftr->Elem1No);
// The test_fe2 object is really a dummy and not used on the
// boundaries, but we can't dereference a NULL pointer, and we don't
fem/bilininteg.cpp
+5 -5
View File
@@ -2949,16 +2949,16 @@ void VectorDiffusionIntegrator::AssembleElementVector(
}
dshape.SetSize(dof, dim);
dshapedxt.SetSize(dof, dim);
// pelmat.SetSize(dim);
dshapedxt.SetSize(dof, sdim);
pelmat.SetSize(dof);
elvect.SetSize(dim*dof);
elvect.SetSize(vdim*dof);
// NOTE: DenseMatrix is in column-major order. This is consistent with
// vectors ordered byNODES. In the resulting DenseMatrix, each column
// corresponds to a particular vdim.
DenseMatrix mat_in(elfun.GetData(), dof, dim);
DenseMatrix mat_out(elvect.GetData(), dof, dim);
DenseMatrix mat_in(elfun.GetData(), dof, vdim);
DenseMatrix mat_out(elvect.GetData(), dof, vdim);
const IntegrationRule *ir = IntRule;
if (ir == NULL)
fem/coefficient.cpp
+2 -2
View File
@@ -1568,7 +1568,7 @@ double ComputeGlobalLpNorm(double p, VectorCoefficient &coeff, ParMesh &pmesh,
#endif
VectorQuadratureFunctionCoefficient::VectorQuadratureFunctionCoefficient(
QuadratureFunction &qf)
const QuadratureFunction &qf)
: VectorCoefficient(qf.GetVDim()), QuadF(qf), index(0) { }
void VectorQuadratureFunctionCoefficient::SetComponent(int index_, int length_)
@@ -1622,7 +1622,7 @@ void VectorQuadratureFunctionCoefficient::Project(QuadratureFunction &qf)
}
QuadratureFunctionCoefficient::QuadratureFunctionCoefficient(
QuadratureFunction &qf) : QuadF(qf)
const QuadratureFunction &qf) : QuadF(qf)
{
MFEM_VERIFY(qf.GetVDim() == 1, "QuadratureFunction's vdim must be 1");
}
fem/coefficient.hpp
+2 -2
View File
@@ -2245,7 +2245,7 @@ private:
public:
/// Constructor with a quadrature function as input
VectorQuadratureFunctionCoefficient(QuadratureFunction &qf);
VectorQuadratureFunctionCoefficient(const QuadratureFunction &qf);
/** Set the starting index within the QuadFunc that'll be used to project
outwards as well as the corresponding length. The projected length should
@@ -2273,7 +2273,7 @@ private:
public:
/// Constructor with a quadrature function as input
QuadratureFunctionCoefficient(QuadratureFunction &qf);
QuadratureFunctionCoefficient(const QuadratureFunction &qf);
const QuadratureFunction& GetQuadFunction() const { return QuadF; }
fem/eltrans.hpp
+7 -4
View File
@@ -170,12 +170,13 @@ public:
virtual int GetSpaceDim() const = 0;
/** @brief Transform a point @a pt from physical space to a point @a ip in
reference space. */
reference space and optionally can set a solver tolerance using @a phys_tol. */
/** Attempt to find the IntegrationPoint that is transformed into the given
point in physical space. If the inversion fails a non-zero value is
returned. This method is not 100 percent reliable for non-linear
transformations. */
virtual int TransformBack(const Vector &pt, IntegrationPoint &ip) = 0;
virtual int TransformBack(const Vector &pt, IntegrationPoint &ip,
const double phys_tol = 1e-15) = 0;
virtual ~ElementTransformation() { }
};
@@ -442,14 +443,16 @@ public:
virtual int GetSpaceDim() const { return PointMat.Height(); }
/** @brief Transform a point @a pt from physical space to a point @a ip in
reference space. */
reference space and optionally can set a solver tolerance using @a phys_tol. */
/** Attempt to find the IntegrationPoint that is transformed into the given
point in physical space. If the inversion fails a non-zero value is
returned. This method is not 100 percent reliable for non-linear
transformations. */
virtual int TransformBack(const Vector & v, IntegrationPoint & ip)
virtual int TransformBack(const Vector & v, IntegrationPoint & ip,
const double phys_rel_tol = 1e-15)
{
InverseElementTransformation inv_tr(this);
inv_tr.SetPhysicalRelTol(phys_rel_tol);
return inv_tr.Transform(v, ip);
}
fem/gridfunc.cpp
+9 -8
View File
@@ -2058,7 +2058,7 @@ void GridFunction::AccumulateAndCountZones(VectorCoefficient &vcoeff,
}
void GridFunction::AccumulateAndCountBdrValues(
Coefficient *coeff[], VectorCoefficient *vcoeff, Array<int> &attr,
Coefficient *coeff[], VectorCoefficient *vcoeff, const Array<int> &attr,
Array<int> &values_counter)
{
int i, j, fdof, d, ind, vdim;
@@ -2200,7 +2200,7 @@ static void accumulate_dofs(const Array<int> &dofs, const Vector &vals,
}
void GridFunction::AccumulateAndCountBdrTangentValues(
VectorCoefficient &vcoeff, Array<int> &bdr_attr,
VectorCoefficient &vcoeff, const Array<int> &bdr_attr,
Array<int> &values_counter)
{
const FiniteElement *fe;
@@ -2576,7 +2576,7 @@ void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
}
void GridFunction::ProjectBdrCoefficient(VectorCoefficient &vcoeff,
Array<int> &attr)
const Array<int> &attr)
{
Array<int> values_counter;
AccumulateAndCountBdrValues(NULL, &vcoeff, attr, values_counter);
@@ -2593,7 +2593,8 @@ void GridFunction::ProjectBdrCoefficient(VectorCoefficient &vcoeff,
#endif
}
void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[], Array<int> &attr)
void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[],
const Array<int> &attr)
{
Array<int> values_counter;
// this->HostReadWrite(); // done inside the next call
@@ -2623,7 +2624,7 @@ void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[], Array<int> &attr)
}
void GridFunction::ProjectBdrCoefficientNormal(
VectorCoefficient &vcoeff, Array<int> &bdr_attr)
VectorCoefficient &vcoeff, const Array<int> &bdr_attr)
{
#if 0
// implementation for the case when the face dofs are integrals of the
@@ -2698,7 +2699,7 @@ void GridFunction::ProjectBdrCoefficientNormal(
}
void GridFunction::ProjectBdrCoefficientTangent(
VectorCoefficient &vcoeff, Array<int> &bdr_attr)
VectorCoefficient &vcoeff, const Array<int> &bdr_attr)
{
Array<int> values_counter;
AccumulateAndCountBdrTangentValues(vcoeff, bdr_attr, values_counter);
@@ -2746,7 +2747,8 @@ double GridFunction::ComputeL2Error(
for (j = 0; j < ir->GetNPoints(); j++)
{
const IntegrationPoint &ip = ir->IntPoint(j);
fe->CalcShape(ip, shape);
transf->SetIntPoint(&ip);
fe->CalcPhysShape(*transf, shape);
for (d = 0; d < fes->GetVDim(); d++)
{
a = 0;
@@ -2759,7 +2761,6 @@ double GridFunction::ComputeL2Error(
{
a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
}
transf->SetIntPoint(&ip);
a -= exsol[d]->Eval(*transf, ip);
error += ip.weight * transf->Weight() * a * a;
}
fem/gridfunc.hpp
+9 -7
View File
@@ -446,11 +446,12 @@ protected:
Array<int> &zones_per_dof);
void AccumulateAndCountBdrValues(Coefficient *coeff[],
VectorCoefficient *vcoeff, Array<int> &attr,
VectorCoefficient *vcoeff,
const Array<int> &attr,
Array<int> &values_counter);
void AccumulateAndCountBdrTangentValues(VectorCoefficient &vcoeff,
Array<int> &bdr_attr,
const Array<int> &bdr_attr,
Array<int> &values_counter);
// Complete the computation of averages; called e.g. after
@@ -465,7 +466,7 @@ public:
/** @brief Project a Coefficient on the GridFunction, modifying only DOFs on
the boundary associated with the boundary attributes marked in the
@a attr array. */
void ProjectBdrCoefficient(Coefficient &coeff, Array<int> &attr)
void ProjectBdrCoefficient(Coefficient &coeff, const Array<int> &attr)
{
Coefficient *coeff_p = &coeff;
ProjectBdrCoefficient(&coeff_p, attr);
@@ -475,26 +476,27 @@ public:
DOFs on the boundary associated with the boundary attributes marked in
the @a attr array. */
virtual void ProjectBdrCoefficient(VectorCoefficient &vcoeff,
Array<int> &attr);
const Array<int> &attr);
/** @brief Project a set of Coefficient%s on the components of the
GridFunction, modifying only DOFs on the boundary associated with the
boundary attributed marked in the @a attr array. */
/** If a Coefficient pointer in the array @a coeff is NULL, that component
will not be touched. */
virtual void ProjectBdrCoefficient(Coefficient *coeff[], Array<int> &attr);
virtual void ProjectBdrCoefficient(Coefficient *coeff[],
const Array<int> &attr);
/** Project the normal component of the given VectorCoefficient on
the boundary. Only boundary attributes that are marked in
'bdr_attr' are projected. Assumes RT-type VectorFE GridFunction. */
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff,
Array<int> &bdr_attr);
const Array<int> &bdr_attr);
/** @brief Project the tangential components of the given VectorCoefficient
on the boundary. Only boundary attributes that are marked in @a bdr_attr
are projected. Assumes ND-type VectorFE GridFunction. */
virtual void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
Array<int> &bdr_attr);
const Array<int> &bdr_attr);
virtual double ComputeL2Error(Coefficient *exsol[],
const IntegrationRule *irs[] = NULL,
fem/intrules.cpp
+1
View File
@@ -917,6 +917,7 @@ int Quadrature1D::CheckClosed(int type)
{
case GaussLobatto:
case ClosedUniform:
case ClosedGL:
return type;
default:
return Invalid;
fem/lininteg.hpp
+2 -2
View File
@@ -728,7 +728,7 @@ private:
public:
VectorQuadratureLFIntegrator(VectorQuadratureFunctionCoefficient &vqfc,
const IntegrationRule *ir)
const IntegrationRule *ir = NULL)
: LinearFormIntegrator(ir), vqfc(vqfc)
{
if (ir)
@@ -760,7 +760,7 @@ private:
public:
QuadratureLFIntegrator(QuadratureFunctionCoefficient &qfc,
const IntegrationRule *ir)
const IntegrationRule *ir = NULL)
: LinearFormIntegrator(ir), qfc(qfc)
{
if (ir)
fem/pgridfunc.cpp
+2 -2
View File
@@ -666,7 +666,7 @@ void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &vcoeff,
}
void ParGridFunction::ProjectBdrCoefficient(
Coefficient *coeff[], VectorCoefficient *vcoeff, Array<int> &attr)
Coefficient *coeff[], VectorCoefficient *vcoeff, const Array<int> &attr)
{
Array<int> values_counter;
AccumulateAndCountBdrValues(coeff, vcoeff, attr, values_counter);
@@ -720,7 +720,7 @@ void ParGridFunction::ProjectBdrCoefficient(
}
void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
Array<int> &bdr_attr)
const Array<int> &bdr_attr)
{
Array<int> values_counter;
AccumulateAndCountBdrTangentValues(vcoeff, bdr_attr, values_counter);
fem/pgridfunc.hpp
+5 -4
View File
@@ -44,7 +44,7 @@ protected:
Vector send_data;
void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff,
Array<int> &attr);
const Array<int> &attr);
public:
ParGridFunction() { pfes = NULL; }
@@ -256,16 +256,17 @@ public:
// Only the values in the master are guaranteed to be correct!
void ProjectBdrCoefficient(VectorCoefficient &vcoeff,
Array<int> &attr) override
const Array<int> &attr) override
{ ProjectBdrCoefficient(NULL, &vcoeff, attr); }
// Only the values in the master are guaranteed to be correct!
void ProjectBdrCoefficient(Coefficient *coeff[], Array<int> &attr) override
void ProjectBdrCoefficient(Coefficient *coeff[],
const Array<int> &attr) override
{ ProjectBdrCoefficient(coeff, NULL, attr); }
// Only the values in the master are guaranteed to be correct!
void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
Array<int> &bdr_attr) override;
const Array<int> &bdr_attr) override;
double ComputeL1Error(Coefficient *exsol[],
const IntegrationRule *irs[] = NULL) const override
fem/transfer.cpp
+42 -12
View File
@@ -1420,14 +1420,20 @@ void Prolongation2D(const int NE, const int D1D, const int Q1D,
const Array<double>& B, const Vector& mask)
{
auto x_ = Reshape(localL.Read(), D1D, D1D, NE);
auto y_ = Reshape(localH.ReadWrite(), Q1D, Q1D, NE);
auto y_ = Reshape(localH.Write(), Q1D, Q1D, NE);
auto B_ = Reshape(B.Read(), Q1D, D1D);
auto m_ = Reshape(mask.Read(), Q1D, Q1D, NE);
localH = 0.0;
mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
{
for (int qy = 0; qy < Q1D; ++qy)
{
for (int qx = 0; qx < Q1D; ++qx)
{
y_(qx, qy, e) = 0.0;
}
}
for (int dy = 0; dy < D1D; ++dy)
{
double sol_x[DofQuadLimits::MAX_Q1D];
@@ -1467,14 +1473,23 @@ void Prolongation3D(const int NE, const int D1D, const int Q1D,
const Array<double>& B, const Vector& mask)
{
auto x_ = Reshape(localL.Read(), D1D, D1D, D1D, NE);
auto y_ = Reshape(localH.ReadWrite(), Q1D, Q1D, Q1D, NE);
auto y_ = Reshape(localH.Write(), Q1D, Q1D, Q1D, NE);
auto B_ = Reshape(B.Read(), Q1D, D1D);
auto m_ = Reshape(mask.Read(), Q1D, Q1D, Q1D, NE);
localH = 0.0;
mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
{
for (int qz = 0; qz < Q1D; ++qz)
{
for (int qy = 0; qy < Q1D; ++qy)
{
for (int qx = 0; qx < Q1D; ++qx)
{
y_(qx, qy, qz, e) = 0.0;
}
}
}
for (int dz = 0; dz < D1D; ++dz)
{
double sol_xy[DofQuadLimits::MAX_Q1D][DofQuadLimits::MAX_Q1D];
@@ -1539,14 +1554,20 @@ void Restriction2D(const int NE, const int D1D, const int Q1D,
const Array<double>& Bt, const Vector& mask)
{
auto x_ = Reshape(localH.Read(), Q1D, Q1D, NE);
auto y_ = Reshape(localL.ReadWrite(), D1D, D1D, NE);
auto y_ = Reshape(localL.Write(), D1D, D1D, NE);
auto Bt_ = Reshape(Bt.Read(), D1D, Q1D);
auto m_ = Reshape(mask.Read(), Q1D, Q1D, NE);
localL = 0.0;
mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
{
for (int dy = 0; dy < D1D; ++dy)
{
for (int dx = 0; dx < D1D; ++dx)
{
y_(dx, dy, e) = 0.0;
}
}
for (int qy = 0; qy < Q1D; ++qy)
{
double sol_x[DofQuadLimits::MAX_D1D];
@@ -1578,14 +1599,23 @@ void Restriction3D(const int NE, const int D1D, const int Q1D,
const Array<double>& Bt, const Vector& mask)
{
auto x_ = Reshape(localH.Read(), Q1D, Q1D, Q1D, NE);
auto y_ = Reshape(localL.ReadWrite(), D1D, D1D, D1D, NE);
auto y_ = Reshape(localL.Write(), D1D, D1D, D1D, NE);
auto Bt_ = Reshape(Bt.Read(), D1D, Q1D);
auto m_ = Reshape(mask.Read(), Q1D, Q1D, Q1D, NE);
localL = 0.0;
mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
{
for (int dz = 0; dz < D1D; ++dz)
{
for (int dy = 0; dy < D1D; ++dy)
{
for (int dx = 0; dx < D1D; ++dx)
{
y_(dx, dy, dz, e) = 0.0;
}
}
}
for (int qz = 0; qz < Q1D; ++qz)
{
double sol_xy[DofQuadLimits::MAX_D1D][DofQuadLimits::MAX_D1D];
linalg/blockoperator.cpp
+2 -2
View File
@@ -45,12 +45,12 @@ BlockOperator::BlockOperator(const Array<int> & row_offsets_,
op = static_cast<Operator *>(NULL);
}
void BlockOperator::SetDiagonalBlock(int iblock, Operator *opt, double c)
void BlockOperator::SetDiagonalBlock(int iblock, const Operator *opt, double c)
{
SetBlock(iblock, iblock, opt, c);
}
void BlockOperator::SetBlock(int iRow, int iCol, Operator *opt, double c)
void BlockOperator::SetBlock(int iRow, int iCol, const Operator *opt, double c)
{
if (owns_blocks && op(iRow, iCol))
{
linalg/blockoperator.hpp
+3 -6
View File
@@ -61,14 +61,14 @@ public:
* op: the Operator to be inserted.
* c: optional scalar multiple for this block.
*/
void SetDiagonalBlock(int iblock, Operator *op, double c = 1.0);
void SetDiagonalBlock(int iblock, const Operator *op, double c = 1.0);
//! Add a block op in the block-entry (iblock, jblock).
/**
* irow, icol: The block will be inserted in location (irow, icol).
* op: the Operator to be inserted.
* c: optional scalar multiple for this block.
*/
void SetBlock(int iRow, int iCol, Operator *op, double c = 1.0);
void SetBlock(int iRow, int iCol, const Operator *op, double c = 1.0);
//! Return the number of row blocks
int NumRowBlocks() const { return nRowBlocks; }
@@ -78,9 +78,6 @@ public:
//! Check if block (i,j) is a zero block
int IsZeroBlock(int i, int j) const { return (op(i,j)==NULL) ? 1 : 0; }
//! Return a reference to block i,j
Operator & GetBlock(int i, int j)
{ MFEM_VERIFY(op(i,j), ""); return *op(i,j); }
//! Return a reference to block i,j (const version)
const Operator & GetBlock(int i, int j) const
{ MFEM_VERIFY(op(i,j), ""); return *op(i,j); }
//! Return the coefficient for block i,j
@@ -123,7 +120,7 @@ private:
//! Column offsets for the starting position of each block
Array<int> col_offsets;
//! 2D array that stores each block of the operator.
Array2D<Operator *> op;
Array2D<const Operator *> op;
//! 2D array that stores a coefficient for each block of the operator.
Array2D<double> coef;
linalg/operator.cpp
+10 -4
View File
@@ -519,7 +519,8 @@ ConstrainedOperator::ConstrainedOperator(Operator *A, const Array<int> &list,
MemoryType mem_type = GetMemoryType(mem_class);
list.Read(); // TODO: just ensure 'list' is registered, no need to copy it
constraint_list.MakeRef(list);
// typically z and w are large vectors, so store them on the device
// typically z and w are large vectors, so use the device (GPU) to perform
// operations on them
z.SetSize(height, mem_type); z.UseDevice(true);
w.SetSize(height, mem_type); w.UseDevice(true);
}
@@ -629,6 +630,13 @@ void ConstrainedOperator::Mult(const Vector &x, Vector &y) const
}
}
void ConstrainedOperator::AddMult(const Vector &x, Vector &y,
const double a) const
{
Mult(x, w);
y.Add(a, w);
}
RectangularConstrainedOperator::RectangularConstrainedOperator(
Operator *A,
const Array<int> &trial_list,
@@ -663,9 +671,7 @@ void RectangularConstrainedOperator::EliminateRHS(const Vector &x,
d_w[id] = d_x[id];
});
// A.AddMult(w, b, -1.0); // if available to all Operators
A->Mult(w, z);
b -= z;
A->AddMult(w, b, -1.0);
const int test_csz = test_constraints.Size();
auto test_idx = test_constraints.Read();
linalg/operator.hpp
+6 -4
View File
@@ -914,14 +914,14 @@ public:
DiagonalPolicy diag_policy = DIAG_ONE);
/// Returns the type of memory in which the solution and temporaries are stored.
virtual MemoryClass GetMemoryClass() const { return mem_class; }
MemoryClass GetMemoryClass() const override { return mem_class; }
/// Set the diagonal policy for the constrained operator.
void SetDiagonalPolicy(const DiagonalPolicy diag_policy_)
{ diag_policy = diag_policy_; }
/// Diagonal of A, modified according to the used DiagonalPolicy.
virtual void AssembleDiagonal(Vector &diag) const;
void AssembleDiagonal(Vector &diag) const override;
/** @brief Eliminate "essential boundary condition" values specified in @a x
from the given right-hand side @a b.
@@ -944,10 +944,12 @@ public:
where the "_b" subscripts denote the essential (boundary) indices/dofs of
the vectors, and "_i" -- the rest of the entries. */
virtual void Mult(const Vector &x, Vector &y) const;
void Mult(const Vector &x, Vector &y) const override;
void AddMult(const Vector &x, Vector &y, const double a = 1.0) const override;
/// Destructor: destroys the unconstrained Operator, if owned.
virtual ~ConstrainedOperator() { if (own_A) { delete A; } }
~ConstrainedOperator() override { if (own_A) { delete A; } }
};
/** @brief Rectangular Operator for imposing essential boundary conditions on
linalg/sparsemat.cpp
+53 -6
View File
@@ -1767,6 +1767,7 @@ void SparseMatrix::EliminateCol(int col, DiagonalPolicy dpolicy)
if (aux->Column == col)
{
aux->Value = 0.0;
break;
}
}
}
@@ -3293,16 +3294,35 @@ void SparseMatrix::PrintMatlab(std::ostream & os) const
{
os << "% size " << height << " " << width << "\n";
os << "% Non Zeros " << NumNonZeroElems() << "\n";
int i, j;
ios::fmtflags old_fmt = os.flags();
os.setf(ios::scientific);
std::streamsize old_prec = os.precision(14);
for (i = 0; i < height; i++)
if (A == NULL)
{
for (j = I[i]; j < I[i+1]; j++)
RowNode *nd;
for (i = 0; i < height; i++)
{
os << i+1 << " " << J[j]+1 << " " << A[j] << '\n';
for (nd = Rows[i], j = 0; nd != NULL; nd = nd->Prev, j++)
{
os << i+1 << " " << nd->Column+1 << " " << nd->Value << '\n';
}
}
}
else
{
// HostRead forces synchronization
HostReadI();
HostReadJ();
HostReadData();
for (i = 0; i < height; i++)
{
for (j = I[i]; j < I[i+1]; j++)
{
os << i+1 << " " << J[j]+1 << " " << A[j] << '\n';
}
}
}
// Write a zero entry at (m,n) to make sure MATLAB doesn't shrink the matrix
@@ -3322,11 +3342,30 @@ void SparseMatrix::PrintMM(std::ostream & os) const
<< "% Generated by MFEM" << '\n';
os << height << " " << width << " " << NumNonZeroElems() << '\n';
for (i = 0; i < height; i++)
if (A == NULL)
{
for (j = I[i]; j < I[i+1]; j++)
RowNode *nd;
for (i = 0; i < height; i++)
{
os << i+1 << " " << J[j]+1 << " " << A[j] << '\n';
for (nd = Rows[i], j = 0; nd != NULL; nd = nd->Prev, j++)
{
os << i+1 << " " << nd->Column+1 << " " << nd->Value << '\n';
}
}
}
else
{
// HostRead forces synchronization
HostReadI();
HostReadJ();
HostReadData();
for (i = 0; i < height; i++)
{
for (j = I[i]; j < I[i+1]; j++)
{
os << i+1 << " " << J[j]+1 << " " << A[j] << '\n';
}
}
}
os.precision(old_prec);
@@ -3341,6 +3380,10 @@ void SparseMatrix::PrintCSR(std::ostream & os) const
os << height << '\n'; // number of rows
// HostRead forces synchronization
HostReadI();
HostReadJ();
HostReadData();
for (i = 0; i <= height; i++)
{
os << I[i]+1 << '\n';
@@ -3366,6 +3409,10 @@ void SparseMatrix::PrintCSR2(std::ostream & os) const
os << height << '\n'; // number of rows
os << width << '\n'; // number of columns
// HostRead forces synchronization
HostReadI();
HostReadJ();
HostReadData();
for (i = 0; i <= height; i++)
{
os << I[i] << '\n';
linalg/sparsemat.hpp
+5
View File
@@ -650,18 +650,23 @@ public:
SparseMatrix &operator*=(double a);
/// Prints matrix to stream out.
/** @note The host in synchronized when the finalized matrix is on the device. */
void Print(std::ostream &out = mfem::out, int width_ = 4) const;
/// Prints matrix in matlab format.
/** @note The host in synchronized when the finalized matrix is on the device. */
virtual void PrintMatlab(std::ostream &out = mfem::out) const;
/// Prints matrix in Matrix Market sparse format.
/** @note The host in synchronized when the finalized matrix is on the device. */
void PrintMM(std::ostream &out = mfem::out) const;
/// Prints matrix to stream out in hypre_CSRMatrix format.
/** @note The host in synchronized when the finalized matrix is on the device. */
void PrintCSR(std::ostream &out) const;
/// Prints a sparse matrix to stream out in CSR format.
/** @note The host in synchronized when the finalized matrix is on the device. */
void PrintCSR2(std::ostream &out) const;
/// Print various sparse matrix statistics.
mesh/pmesh.hpp
+4
View File
@@ -491,6 +491,10 @@ public:
void GenerateOffsets(int N, HYPRE_BigInt loc_sizes[],
Array<HYPRE_BigInt> *offsets[]) const;
/** Return true if the face is interior or shared. In parallel, this
method only works if the face neighbor data is exchanged. */
inline bool FaceIsTrueInterior(int FaceNo) const { return Mesh::FaceIsTrueInterior(FaceNo); }
void ExchangeFaceNbrData();
void ExchangeFaceNbrNodes();
miniapps/solvers/div_free_solver.cpp
+35 -3
View File
@@ -311,6 +311,38 @@ void SaddleSchwarzSmoother::Mult(const Vector & x, Vector & y) const
blk_y.GetBlock(1) -= coarse_l2_projection;
}
BDPMinresSolver::BDPMinresSolver(const HypreParMatrix& M,
const HypreParMatrix& B,
IterSolveParameters param)
: DarcySolver(M.NumRows(), B.NumRows()), op_(offsets_), prec_(offsets_),
BT_(B.Transpose()), solver_(M.GetComm())
{
op_.SetBlock(0,0, &M);
op_.SetBlock(0,1, BT_.As<HypreParMatrix>());
op_.SetBlock(1,0, &B);
Vector Md;
M.GetDiag(Md);
BT_.As<HypreParMatrix>()->InvScaleRows(Md);
S_.Reset(ParMult(&B, BT_.As<HypreParMatrix>()));
BT_.As<HypreParMatrix>()->ScaleRows(Md);
prec_.SetDiagonalBlock(0, new HypreDiagScale(M));
prec_.SetDiagonalBlock(1, new HypreBoomerAMG(*S_.As<HypreParMatrix>()));
static_cast<HypreBoomerAMG&>(prec_.GetDiagonalBlock(1)).SetPrintLevel(0);
prec_.owns_blocks = true;
SetOptions(solver_, param);
solver_.SetOperator(op_);
solver_.SetPreconditioner(prec_);
}
void BDPMinresSolver::Mult(const Vector & x, Vector & y) const
{
solver_.Mult(x, y);
for (int dof : ess_zero_dofs_) { y[dof] = 0.0; }
}
DivFreeSolver::DivFreeSolver(const HypreParMatrix &M, const HypreParMatrix& B,
const DFSData& data)
: DarcySolver(M.NumRows(), B.NumRows()), data_(data), param_(data.param),
@@ -326,8 +358,8 @@ DivFreeSolver::DivFreeSolver(const HypreParMatrix &M, const HypreParMatrix& B,
for (int l = data.P_l2.size(); l >= 0; --l)
{
auto& M_f = static_cast<HypreParMatrix&>(ops_[l]->GetBlock(0, 0));
auto& B_f = static_cast<HypreParMatrix&>(ops_[l]->GetBlock(1, 0));
auto& M_f = static_cast<const HypreParMatrix&>(ops_[l]->GetBlock(0, 0));
auto& B_f = static_cast<const HypreParMatrix&>(ops_[l]->GetBlock(1, 0));
if (l == 0)
{
@@ -552,7 +584,7 @@ void DivFreeSolver::Mult(const Vector & x, Vector & y) const
ch.Clear();
ch.Start();
auto M = dynamic_cast<HypreParMatrix&>(ops_.Last()->GetBlock(0, 0));
auto& M = dynamic_cast<const HypreParMatrix&>(ops_.Last()->GetBlock(0, 0));
M.Mult(-1.0, correction.GetBlock(0), 1.0, resid.GetBlock(0));
SolvePotential(resid.GetBlock(0), correction.GetBlock(1));
blk_y.GetBlock(1) += correction.GetBlock(1);
miniapps/solvers/div_free_solver.hpp
+18
View File
@@ -158,6 +158,24 @@ public:
virtual void SetOperator(const Operator &op) { }
};
/// Wrapper for the block-diagonal-preconditioned MINRES defined in ex5p.cpp
class BDPMinresSolver : public DarcySolver
{
BlockOperator op_;
BlockDiagonalPreconditioner prec_;
OperatorPtr BT_;
OperatorPtr S_; // S_ = B diag(M)^{-1} B^T
MINRESSolver solver_;
Array<int> ess_zero_dofs_;
public:
BDPMinresSolver(const HypreParMatrix& M, const HypreParMatrix& B,
IterSolveParameters param);
virtual void Mult(const Vector & x, Vector & y) const;
virtual void SetOperator(const Operator &op) { }
void SetEssZeroDofs(const Array<int>& dofs) { dofs.Copy(ess_zero_dofs_); }
virtual int GetNumIterations() const { return solver_.GetNumIterations(); }
};
/// Divergence free solver.
/** Divergence free solver.
The basic idea of the solver is to exploit a multilevel decomposition of
tests/unit/linalg/test_matrix_sparse.cpp
+163
View File
@@ -11,6 +11,7 @@
#include "mfem.hpp"
#include "unit_tests.hpp"
#include <sstream>
namespace mfem
{
@@ -77,4 +78,166 @@ TEST_CASE("SparseMatrixAbsMult", "[SparseMatrixAbsMult]")
}
}
TEST_CASE("SparseMatrix printing", "[SparseMatrix]")
{
// Create a test sparse matrix and print it using different methods
// and compare the output with the reference one
DenseMatrix dense(
{
{0.0, 4.0, 0.0},
{5.0, 0.0, 1.0},
{2.0, 0.0, 0.0}
});
const int width = dense.Width();
const int height = dense.Height();
int nonzero = 0;
// Non-finalized matrix (LIL)
SparseMatrix mat_lil(height, width);
for (int i = 0; i < height; i++)
for (int j = 0; j < width; j++)
{
if (dense(i,j) != 0.0)
{
mat_lil.Add(i, j, dense(i,j));
nonzero++;
}
}
// Finalized matrix (CSR)
SparseMatrix mat_csr(mat_lil);
mat_csr.Finalize();
std::stringstream ss, ss_ref;
SECTION("Print")
{
//assume print width >= matrix width
ss_ref.str("");
for (int i = 0; i < height; i++)
{
ss_ref << "[row " << i << "]\n";
for (int j = width-1; j >=0 ; j--)
if (dense(i,j) != 0.0)
{
ss_ref << " (" << j << "," << dense(i,j) << ")";
}
ss_ref << "\n";
}
ss.str("");
mat_lil.Print(ss);
REQUIRE(ss.str() == ss_ref.str());
ss.str("");
mat_csr.Print(ss);
REQUIRE(ss.str() == ss_ref.str());
}
SECTION("PrintMatlab")
{
ss_ref.str("");
ss_ref << "% size " << height << " " << width << "\n";
ss_ref << "% Non Zeros " << nonzero << "\n";
std::ios::fmtflags old_fmt = ss_ref.flags();
ss_ref.setf(std::ios::scientific);
std::streamsize old_prec = ss_ref.precision(14);
for (int i = 0; i < height; i++)
for (int j = width-1; j >=0 ; j--)
if (dense(i,j) != 0.0)
{
ss_ref << i+1 << " " << j+1 << " " << dense(i,j) << "\n";
}
ss_ref << height << " " << width << " 0.0\n";
ss_ref.precision(old_prec);
ss_ref.flags(old_fmt);
ss.str("");
mat_lil.PrintMatlab(ss);
REQUIRE(ss.str() == ss_ref.str());
ss.str("");
mat_csr.PrintMatlab(ss);
REQUIRE(ss.str() == ss_ref.str());
}
SECTION("PrintMM")
{
ss_ref.str("");
ss_ref << "%%MatrixMarket matrix coordinate real general" << '\n'
<< "% Generated by MFEM" << '\n';
ss_ref << height << " " << width << " " << nonzero << "\n";
std::ios::fmtflags old_fmt = ss_ref.flags();
ss_ref.setf(std::ios::scientific);
std::streamsize old_prec = ss_ref.precision(14);
for (int i = 0; i < height; i++)
for (int j = width-1; j >=0 ; j--)
if (dense(i,j) != 0.0)
{
ss_ref << i+1 << " " << j+1 << " " << dense(i,j) << "\n";
}
ss_ref.precision(old_prec);
ss_ref.flags(old_fmt);
ss.str("");
mat_lil.PrintMM(ss);
REQUIRE(ss.str() == ss_ref.str());
ss.str("");
mat_csr.PrintMM(ss);
REQUIRE(ss.str() == ss_ref.str());
}
SECTION("PrintCSR")
{
ss_ref.str("");
ss_ref << height << "\n";
Array<int> I(height+1);
Array<int> J(nonzero);
Vector A(nonzero);
int idx = 0;
for (int i = 0; i < height; i++)
{
I[i] = idx;
for (int j = width-1; j >=0 ; j--)
if (dense(i,j) != 0.0)
{
J[idx] = j;
A[idx] = dense(i,j);
idx++;
}
}
I[height] = idx;
for (int i = 0; i <= height; i++)
{
ss_ref << I[i]+1 << '\n';
}
for (int i = 0; i < I[height]; i++)
{
ss_ref << J[i]+1 << '\n';
}
for (int i = 0; i < I[height]; i++)
{
ss_ref << A[i] << '\n';
}
ss.str("");
mat_csr.PrintCSR(ss);
REQUIRE(ss.str() == ss_ref.str());
}
}
} // namespace mfem