[petsc-users] Construct Matrix based on row and column values
Elias Karabelas
karabelaselias at gmail.com
Wed Apr 8 12:06:04 CDT 2020
Happy to, but this is far from a production code speaking of commenting
and stuff :D
##############################################################################
void CalculateFluxCorrection(ConvectionDiffusionOperator * convOp)
{
//convOp->R SparseMatrix
//convOp->MC consistent Mass Matrix
//convOp->L Fluxcorrected stiffnessmatrix
//convOp->D Helper matrix
//convOp->ML lumped Mass matrix (Vec)
//Zero Matrix R
MatZeroEntries(&convOp->R);
Vec Rpos, Rneg;
VecDuplicate(convOp->utilde, &Rpos);
VecDuplicate(convOp->utilde, &Rneg);
// for later usage we need the maxima and minima of utilde
double utildemax = 0.;
double utildemin = 0.;
VecMax(convOp->utilde, nullptr, &utildemax);
VecMin(convOp->utilde, nullptr, &utildemin);
//Create Scatter to all and globalize the vectors
Vec _vseq, _useq;
VecScatter toall;
VecScatterCreateToAll(convOp->utilde,&toall,&_useq);
VecDuplicate(_useq, &_vseq);
VecScatterBegin(toall, convOp->utilde, _useq, INSERT_VALUES,
SCATTER_FORWARD);
VecScatterEnd(toall, convOp->utilde, _useq, INSERT_VALUES,
SCATTER_FORWARD);
VecScatterBegin(toall, convOp->v, _vseq, INSERT_VALUES, SCATTER_FORWARD);
VecScatterEnd(toall, convOp->v, _vseq, INSERT_VALUES, SCATTER_FORWARD);
//Values of utilde and v globalized on each rank
const double *_useq_ptr, *_vseq_ptr;
VecGetArrayRead(_vseq, &_vseq_ptr);
VecGetArrayRead(_useq, &_useq_ptr);
double *lmassloc;
VecGetArrayRead(convOp->ML, &lmassloc);
//Need write access to local part of Rpos and Rneg since I update the
values
double *rposloc;
VecGetArray(Rpos, &rposloc);
double *rnegloc;
VecGetArray(Rneg, &rnegloc);
//Stuff needed for assembling convOp->R
int rstart, rend;
int ncolsM, ncolsD, ncolsR;
const int *colsM, *colsD, *colsR;
const double *valsM, *valsD, *valsR;
MatGetOwnershipRange(convOp->R,&rstart,&rend);
const size_t sz = rend - rstart;
std::vector<double> Lvals;
//R is defined as R_ij = dt * convOp->MC[i,j] * (v[i] - v[j]) -
convOp->D[i,j] * (utilde[i] - utilde[j])
for(size_t i=0; i < sz; i++)
{
int row = rstart + i;
int diag = 0;
//Fill i-th row of convOp->R
MatGetRow(convOp->MC, row, &ncolsM, &colsM, &valsM);
MatGetRow(convOp->D, row, &ncolsD, &colsD, &valsD);
assert(ncolsM == ncolsD); //Assert symmetric structure should be
given for FEM matrices
Lvals.assign(ncolsM, 0.0);
double Pipos = 0.;
double Pineg = 0.;
for(int c=0; c < ncolsM; c++)
{
assert(colsM[c] == colsD[c]);
if(colsM[c] == row)
diag = c;
else
{
Lvals[c] = convOp->dt * (valsM[c] * (_vseq_ptr[row] -
_vseq_ptr[colsM[c]]) - valsD[c] * (_useq_ptr[row] - _useq_ptr[colsM[c]]));
Pipos += std::max(0.0, Lvals[c]);
Pineg += std::min(0.0, Lvals[c]);
}
}
Lvals[diag] = 0.0;
MatSetValues(convOp->R, 1, &row, ncolsM, colsM, Lvals.data(),
ADD_VALUES);
//Qip[i] := max(0, max_{j=0,...,N-1, j neq i}(utilde[j] - utilde[i]))
//Qim[i] := min(0, min_{j=0,...,N-1, j neq i}(utilde[j] - utilde[i]))
const double Qip = std::max(0.0, utildemax - _useq_ptr[row]);
const double Qim = std::min(0.0, utildemin - _useq_ptr[row]);
//calc R+ and R-
rposloc[i] = Pipos != 0 ? std::min(1.0, lmassloc[i] * Qip / Pipos)
: 0.0;
rnegloc[i] = Pineg != 0 ? std::min(1.0, lmassloc[i] * Qim / Pineg)
: 0.0;
MatRestoreRow(convOp->MC.m, row, &ncolsM, &colsM, &valsM);
MatRestoreRow(convOp->D.m, row, &ncolsD, &colsD, &valsD);
}
MatAssemblyBegin(convOp->R.m,MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(convOp->R.m,MAT_FINAL_ASSEMBLY);
VecRestoreArray(Rpos, &rposloc);
VecRestoreArray(Rneg, &rnegloc);
VecRestoreArrayRead(convOp->ML, &lmassloc);
VecRestoreArrayRead(_vseq, &_vseq_ptr);
VecRestoreArrayRead(_useq, &_useq_ptr);
// Done assembling R
// Calculate flux correction defined via fstar[i] := Sum_j alpha[i,j]
* convOp->R[i,j]
// with alpha[i,j] = R[i,j] > 0 ? min(R+[i], R-[j]) : min(R-[i], Rp[j])
Vec _rpseq, _rmseq;
VecDuplicate(_useq, &_rpseq);
VecDuplicate(_useq, &_rmseq);
VecDestroy(&_useq);
VecDestroy(&_vseq);
VecScatterBegin(toall, Rpos, _rpseq, INSERT_VALUES, SCATTER_FORWARD);
VecScatterEnd(toall, Rpos, _rpseq, INSERT_VALUES, SCATTER_FORWARD);
VecScatterBegin(toall, Rneg, _rmseq, INSERT_VALUES, SCATTER_FORWARD);
VecScatterEnd(toall, Rneg, _rmseq, INSERT_VALUES, SCATTER_FORWARD);
const double *_rpseq_ptr, *_rmseq_ptr;
VecGetArrayRead(_rpseq, &_rpseq_ptr);
VecGetArrayRead(_rmseq, &_rmseq_ptr);
double *fluxcorr;
VecGetArray(convOp->fstar, &fluxcorr);
double aij = 0.;
//Calculate the flux correction
for(size_t i=0; i < sz; i++)
{
int row = rstart + i;
MatGetRow(convOp->R.m, row, &ncolsR, &colsR, &valsR);
for(int c=0; c < ncolsR; c++)
{
if(colsR[c] != row)
{
aij = valsR[c] > 0 ? std::min(_rpseq_ptr[row],
_rmseq_ptr[colsR[c]]) : std::min(_rmseq_ptr[row], _rpseq_ptr[colsR[c]]);
fluxcorr[i] += aij * valsR[c];
}
}
MatRestoreRow(convOp->R.m, row, &ncolsR, &colsR, &valsR);
}
VecRestoreArray(convOp->fstarm &fluxcorr);
VecRestoreArrayRead(_rpseq, &_rpseq_ptr);
VecRestoreArrayRead(_rmseq, &_rmseq_ptr);
VecDestroy(&_rpseq);
VecDestroy(&_rmseq);
//Cleanup
VecDestroy(Rpos);
VecDestroy(Rneg);
VecScatterDestroy(&toall);
}
On 08/04/2020 18:26, Jed Brown wrote:
> Elias Karabelas <karabelaselias at gmail.com> writes:
>
>> Dear Jed,
>>
>> I'm done implementing my FCT-Solver and it works fairly well on a small
>> amount of MPI-Procs. Additionally to your little snippet I have used a
>> VecScatterToAll.
> Can you share the code including that VecScatterToAll?
>
> There more local ways to get that information.
>
>> Reason is that the flux correction f_i takes the form Sum_{j} alpha_ij
>> r_ij where r_ij is defined on the sparsity pattern of my FEM Matrix and
>> alpha_ij is based on two vectors Rp and Rm. So basically I need
>> off-process values of these vectors to construct alpha which I made with
>> a VecScatterToAll. However I guess this will slow down my overall
>> program quite significant. Any Ideas?
>>
>> Best regards
>>
>> Elias
>>
>> On 23/03/2020 15:53, Jed Brown wrote:
>>> Thanks; please don't drop the list.
>>>
>>> I'd be curious whether this operation is common enough that we should
>>> add it to PETSc. My hesitance has been that people may want many
>>> different variants when working with systems of equations, for example.
>>>
>>> Elias Karabelas <karabelaselias at gmail.com> writes:
>>>
>>>> Dear Jed,
>>>>
>>>> Yes the Matrix A comes from assembling a FEM-convection-diffusion
>>>> operator over a tetrahedral mesh. So my matrix graph should be
>>>> symmetric. Thanks for the snippet
>>>>
>>>> On 23/03/2020 15:42, Jed Brown wrote:
>>>>> Elias Karabelas <karabelaselias at gmail.com> writes:
>>>>>
>>>>>> Dear Users,
>>>>>>
>>>>>> I want to implement a FCT (flux corrected transport) scheme with PETSc.
>>>>>> To this end I have amongst other things create a Matrix whose entries
>>>>>> are given by
>>>>>>
>>>>>> L_ij = -max(0, A_ij, A_ji) for i neq j
>>>>>>
>>>>>> L_ii = Sum_{j=0,..n, j neq i} L_ij
>>>>>>
>>>>>> where Mat A is an (non-symmetric) Input Matrix created beforehand.
>>>>>>
>>>>>> I was wondering how to do this. My first search brought me to
>>>>>> https://www.mcs.anl.gov/petsc/petsc-current/src/mat/examples/tutorials/ex16.c.html
>>>>>>
>>>>>>
>>>>>> but this just goes over the rows of one matrix to set new values and now
>>>>>> I would need to run over the rows and columns of the matrix. My Idea was
>>>>>> to just create a transpose of A and do the same but then the row-layout
>>>>>> will be different and I can't use the same for loop for A and AT and
>>>>>> thus also won't be able to calculate the max's above.
>>>>> Does your matrix have symmetric nonzero structure? (It's typical for
>>>>> finite element methods.)
>>>>>
>>>>> If so, all the indices will match up so I think you can do something like:
>>>>>
>>>>> for (row=rowstart; row<rowend; row++) {
>>>>> PetscScalar Lvals[MAX_LEN];
>>>>> PetscInt diag;
>>>>> MatGetRow(A, row, &ncols, &cols, &vals);
>>>>> MatGetRow(At, row, &ncolst, &colst, &valst);
>>>>> assert(ncols == ncolst); // symmetric structure
>>>>> PetscScalar sum = 0;
>>>>> for (c=0; c<ncols; c++) {
>>>>> assert(cols[c] == colst[c]); // symmetric structure
>>>>> if (cols[c] == row) diag = c;
>>>>> else sum -= (Lvals[c] = -max(0, vals[c], valst[c]));
>>>>> }
>>>>> Lvals[diag] = sum;
>>>>> MatSetValues(L, 1, &row, ncols, cols, Lvals, INSERT_VALUES);
>>>>> MatRestoreRow(A, row, &ncols, &cols, &vals);
>>>>> MatRestoreRow(At, row, &ncolst, &colst, &valst);
>>>>> }
More information about the petsc-users
mailing list