hi,I've been trying creating a matrix with DMCreateMatrix() and then adding extra blocks of nonzeros into it using MatSetValuesBlocked(), but getting some unexpected results if I set the matrix type to BAIJ. It seems to behave as expected if I use matrix type AIJ.
I've attached a minimal example program. It reads in the DMPlex from file, sets up a section on it, creates a matrix (blocksize 2) and then inserts a single 2x2 block at global block indices (0,7). It views the matrix before and after the insertion.
If I run with "-dm_mat_type aij" it gives the expected results, but with "-dm_mat_type baij" it doesn't - e.g. if run in serial, it adds the new nonzeros in the right place but also adds a whole lot of other duplicated entries in block row 0.
Is there something I'm not understanding about BAIJ, or about MatSetValuesBlocked()? or possibly some other mistake?
- Adrian On 20/05/24 12:24 pm, Barry Smith wrote:
You can call MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR) then insert the new values. If it is just a handful of new insertions the extra time should be small.Making a copy of the matrix won't give you a new matrix that is any faster to insert into so best to just use the same matrix.BarryOn May 19, 2024, at 7:44 PM, Adrian Croucher <[email protected]> wrote:This Message Is From an External Sender This message came from outside your organization. hi, I have a Jacobian matrix created using DMCreateMatrix(). What would be the best way to add extra nonzero entries into it? I'm guessing that DMCreateMatrix() allocates the storage so the nonzero structure can't really be easily modified. Would it be a case of creating a new matrix, copying the nonzero entries from the original one and then adding the extra ones, before calling MatSetUp() or similar? If so, how exactly would you copy the nonzero structure from the original matrix? Background: the flow problem I'm solving (on a DMPlex with finite volume method) has complex source terms that depend on the solution (e.g. pressure), and can also depend on other source terms. A simple example is when fluid is extracted from one location, with a pressure-dependent flow rate, and some of it is then reinjected in another location. This can result in poor nonlinear solver convergence. I think the reason is that there are effectively missing Jacobian entries in the row for the reinjection cell, which should have an additional dependence on the solution in the cell where fluid is extracted. - Adrian -- Dr Adrian Croucher Senior Research Fellow Department of Engineering Science Waipapa Taumata Rau / University of Auckland, New Zealand email:[email protected] tel: +64 (0)9 923 4611
-- Dr Adrian Croucher Senior Research Fellow Department of Engineering Science Waipapa Taumata Rau / University of Auckland, New Zealand email:[email protected] tel: +64 (0)9 923 4611
3x3grid.exo
Description: Binary data
program matmodify
#include <petsc/finclude/petsc.h>
use petsc
implicit none
PetscInt, parameter :: overlap = 1
DM :: dm, dist_dm
PetscErrorCode :: ierr
PetscSection :: section
PetscReal, parameter :: vals(4) = [-1.d0, -2.d0, -3.d0, -4.d0]
Mat :: J
call PetscInitialize(PETSC_NULL_CHARACTER, ierr); CHKERRQ(ierr)
call DMPlexCreateFromFile(PETSC_COMM_WORLD, "3x3grid.exo", "mesh", &
PETSC_TRUE, dm, ierr); CHKERRQ(ierr)
call DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE, ierr)
call DMSetFromOptions(dm, ierr)
call DMPlexDistribute(dm, overlap, PETSC_NULL_SF, dist_dm, ierr)
CHKERRQ(ierr)
if (dist_dm .ne. PETSC_NULL_DM) then
call DMDestroy(dm, ierr); CHKERRQ(ierr)
dm = dist_dm
end if
call dm_set_fields(dm, num_components = [1,1])
section = dm_create_section(dm, num_components = [1,1], field_dim = [3,3])
call DMSetSection(dm, section, ierr); CHKERRQ(ierr)
call DMCreateMatrix(dm, J, ierr)
call MatView(J, PETSC_VIEWER_STDOUT_WORLD, ierr)
call MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE, ierr)
call MatSetValuesBlocked(J, 1, [0], 1, [7], vals, INSERT_VALUES, ierr)
call MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY, ierr); CHKERRQ(ierr)
call MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY, ierr); CHKERRQ(ierr)
call MatView(J, PETSC_VIEWER_STDOUT_WORLD, ierr)
call MatDestroy(J, ierr)
call DMDestroy(dm, ierr); CHKERRQ(ierr)
call PetscFinalize(ierr); CHKERRQ(ierr)
contains
subroutine dm_set_fields(dm, num_components)
!! Sets fields in the DM.
DM, intent(in out) :: dm
PetscInt, intent(in) :: num_components(:) !! Number of components in each field
! Locals:
PetscInt :: dim, f
PetscFV :: fvm
PetscErrorCode :: ierr
call DMGetDimension(dm, dim, ierr); CHKERRQ(ierr)
call DMClearFields(dm, ierr); CHKERRQ(ierr)
associate(num_fields => size(num_components))
do f = 1, num_fields
call PetscFVCreate(PETSC_COMM_WORLD, fvm, ierr); CHKERRQ(ierr)
call PetscFVSetFromOptions(fvm, ierr); CHKERRQ(ierr)
call PetscFVSetNumComponents(fvm, num_components(f), ierr); CHKERRQ(ierr)
call PetscFVSetSpatialDimension(fvm, dim, ierr); CHKERRQ(ierr)
call DMAddField(dm, PETSC_NULL_DMLABEL, fvm, ierr); CHKERRQ(ierr)
call PetscFVDestroy(fvm, ierr); CHKERRQ(ierr)
end do
end associate
call DMCreateDS(dm, ierr); CHKERRQ(ierr)
end subroutine dm_set_fields
PetscSection function dm_create_section(dm, num_components, field_dim) result(section)
!! Creates section from the given DM and data layout parameters.
DM, intent(in) :: dm !! DM object
PetscInt, target, intent(in) :: num_components(:) !! Number of components in each field
PetscInt, intent(in) :: field_dim(:) !! Dimension each field is defined on (0 = nodes, etc.)
! Locals:
PetscInt :: dim
PetscInt :: i, num_bc
PetscInt, allocatable, target :: num_dof(:)
PetscInt, target :: bc_field(1)
IS, target :: bc_comps(1), bc_points(1)
PetscInt, pointer :: pnum_components(:), pnum_dof(:), pbc_field(:)
IS, pointer :: pbc_comps(:), pbc_points(:)
DMLabel, pointer :: plabel(:)
PetscErrorCode :: ierr
call DMGetDimension(dm, dim, ierr); CHKERRQ(ierr)
associate (num_fields => size(num_components))
allocate(num_dof(num_fields*(dim+1)))
num_dof = 0
do i = 1, num_fields
num_dof((i-1) * (dim+1) + field_dim(i) + 1) = num_components(i)
end do
! Boundary conditions (none):
num_bc = 0
bc_field(1) = 0
pnum_components => num_components
pnum_dof => num_dof
pbc_field => bc_field
pbc_comps => bc_comps
pbc_points => bc_points
plabel => NULL()
call DMPlexCreateSection(dm, plabel, pnum_components, &
pnum_dof, num_bc, pbc_field, pbc_comps, pbc_points, &
PETSC_NULL_IS, section, ierr); CHKERRQ(ierr)
deallocate(num_dof)
end associate
end function dm_create_section
end program matmodify
