Dear Jed,
Thanks I will try to keep you updated. So I managed to get my stuff
running with one process but it fails with more than one. The matrix
assembly as you described it works fine but I forgot that in the FCT
algorithm one also needs to assemble a matrix D that looks kinda like
D_ij = L_ij * (u[i] - u[j])
where u is a vector with solution data from a previous time step. So the
problem is that I have to now access off-process values of the Vector u,
and doing this for each row of D. My basic first idea was to utilize the
VecScatter Example from the Petsc Manual to get the off-process values
per row of L to build this matrix D. But somehow this does not really
work. Any ideas on how to scatter the Vector-data so that I can get the
correct Vector values ?
At the moment I'm using something along this lines (in a loop where I
get the local array of u with VecGetArray --> u_arr)
MatGetRow(L, row, &ncolsL, &colsL, &valsL);
double * Dvals = new double[ncolsM]();
int * idx_to = new int[ncolsL]();
for(int k=0; k < ncolsM; k++)
idx_to[k] = k;
MPI_Barrier(PETSC_COMM_WORLD);
//Scatter the right stuff
Vec x;
VecScatter scatter;
IS from, to;
double *vals;
VecCreateSeq(PETSC_COMM_SELF, ncolsL, &x);
ISCreateGeneral(PETSC_COMM_SELF,ncolsL,colsL,PETSC_COPY_VALUES,&from);
ISCreateGeneral(PETSC_COMM_SELF,ncolsL,idx_to,PETSC_COPY_VALUES,&to);
VecScatterCreate(u,from,x,to,&scatter);
VecScatterBegin(scatter,u,x,INSERT_VALUES,SCATTER_FORWARD);
VecScatterEnd(scatter,u,x,INSERT_VALUES,SCATTER_FORWARD);
VecGetArray(x,&vals);
for(int c=0; c < ncolsL; c++)
{
Dvals[c] = valsL[c] * (vals[c] - u_arr[row])
}
MatSetValues(D, 1, &row, ncolsL, colsL, Dvals, ADD_VALUES);
and so on...
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 <karabelasel...@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 <karabelasel...@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);
}
