BsrMatrix usage error in MPM solver

[Cucchi01]

I think the usage of BsrMatrix.nnz is correct in solve_rheology.py and solver_implicit_mpm.py .
In particular, BsrMatrix are passed to solve_rheology and then used for:

    # Collider contacts
    subgrid_collisions = collider_mat.nnz > 0
    if subgrid_collisions:
    	...

and

    if compliance_mat.nnz == 0:
        compliance_mat_diagonal = None
    else:
        compliance_mat_diagonal = _register_temp(fem.borrow_temporary(temporary_store, shape=stress.shape, dtype=mat66))
        ....

If I am not wrong, the only operations that sets them are inside fem.integrate or fem.interpolate, and they eventually get defined by bsr_set_from_triplets, e.g.:

    if output is None or add_to_output:
        bsr_result = bsr_zeros(
            rows_of_blocks=test.space_partition.node_count(),
            cols_of_blocks=trial.space_partition.node_count(),
            block_type=block_type,
            device=device,
        )
    else:
        bsr_result = output

    bsr_set_from_triplets(bsr_result, triplet_rows, triplet_cols, triplet_values, **(bsr_options or {}))

In Warp sparse.py, it seems that depending if nnz_sync is called or not, nnz can have the role of capacity or the actual 'number of non-zero'.
In the case of bsr_set_from_triplets, this sync is not called (I think to save synchronization time) and nnz represents the capacity. While, when .nnz is called in the MPM, semantically it seems to need 'number of non-zero'.

A fix could be use the last cell of offsets instead.

Maybe for the API of BsrMatrix an alternative design could be to rename the current nnz and related method to 'capacity', and to include nnz as a @Property that represents the actual 'number of non-zero'.

Please correct if I am wrong, thank you in advance for the help!

[gdaviet]

Hi @Cucchi01,

The nnz field in warp's BsrMatrix is not really a capacity; it is rather an upper bound for the maximum possible number of non-zeros as seen by the host. When doing bsr_set_from_triplets, the host knows that there cannot be more non-zeros than the number of triplets (though there could be less if some triplets are repeated, as indicated by the last element of the device offsets array). When doing bsr_set_zero, the host knows that the nnz sill be exactly zero, so the nnz field will reflect that. On the other hand, the capacity of the matrix is not changed (the arrays are not re-allocated).

As far as the MPM solver in Newton is concerned, the compliance/subgrid_collisions matrices are either set to zero or set as the result of an integration routine, depending on solver parameters. In the first case, the nnz field will reflect that, do the check looks correct to me. For the second case, the nnz>0 branch should still support empty matrices anyway.

That being said, I remember there being a bug where subgrid collisions were always active; I'll check if that has been fixed in the main branch, or just in my local one