1212from devito .passes .iet .languages .C import CDataManager
1313from devito .petsc .types import (DM , Mat , Vec , PetscMPIInt , KSP ,
1414 PC , KSPConvergedReason , PETScArray ,
15- LinearSolveExpr , FieldData , MultipleFieldData )
15+ LinearSolveExpr , FieldData , MultipleFieldData ,
16+ SubMatrixBlock )
1617from devito .petsc .solve import PETScSolve , EssentialBC
1718from devito .petsc .iet .nodes import Expression
1819from devito .petsc .initialize import PetscInitialize
@@ -647,7 +648,7 @@ def define(self, dimensions):
647648 assert jac .col_target == e
648649
649650 # 2 symbolic expressions for each each EssentialBC (One ZeroRow and one ZeroColumn).
650- # NOTE: this is likely to change when PetscSection + DMDA is supported
651+ # NOTE: This is likely to change when PetscSection + DMDA is supported
651652 assert len (jac .matvecs ) == 5
652653 # TODO: I think some internals are preventing symplification here?
653654 assert str (jac .scdiag ) == 'h_x*(1 - 2.0/h_x**2)'
@@ -656,6 +657,48 @@ def define(self, dimensions):
656657 assert not isinstance (jac .matvecs [- 1 ], EssentialBC )
657658
658659
660+ @skipif ('petsc' )
661+ def test_residual ():
662+ class SubLeft (SubDomain ):
663+ name = 'subleft'
664+
665+ def define (self , dimensions ):
666+ x , = dimensions
667+ return {x : ('left' , 1 )}
668+
669+ class SubRight (SubDomain ):
670+ name = 'subright'
671+
672+ def define (self , dimensions ):
673+ x , = dimensions
674+ return {x : ('right' , 1 )}
675+
676+ sub1 = SubLeft ()
677+ sub2 = SubRight ()
678+
679+ grid = Grid (shape = (11 ,), subdomains = (sub1 , sub2 ), dtype = np .float64 )
680+
681+ e = Function (name = 'e' , grid = grid , space_order = 2 )
682+ f = Function (name = 'f' , grid = grid , space_order = 2 )
683+
684+ bc_1 = EssentialBC (e , 1.0 , subdomain = sub1 )
685+ bc_2 = EssentialBC (e , 2.0 , subdomain = sub2 )
686+
687+ eq1 = Eq (e .laplace + e , f + 2.0 )
688+
689+ petsc = PETScSolve ([eq1 , bc_1 , bc_2 ], target = e )
690+
691+ res = petsc .rhs .fielddata .residual
692+
693+ assert res .target == e
694+ # NOTE: This is likely to change when PetscSection + DMDA is supported
695+ assert len (res .F_exprs ) == 5
696+ assert len (res .b_exprs ) == 3
697+
698+ assert not res .time_mapper
699+ assert str (res .scdiag ) == 'h_x*(1 - 2.0/h_x**2)'
700+
701+
659702class TestCoupledLinear :
660703 # The coupled interface can be used even for uncoupled problems, meaning
661704 # the equations will be solved within a single matrix system.
@@ -872,6 +915,42 @@ def test_mixed_jacobian(self):
872915 j10 = jacobian .get_submatrix (1 , 0 )
873916 j11 = jacobian .get_submatrix (1 , 1 )
874917
918+ # Check type of each submatrix is a SubMatrixBlock
919+ assert isinstance (j00 , SubMatrixBlock )
920+ assert isinstance (j01 , SubMatrixBlock )
921+ assert isinstance (j10 , SubMatrixBlock )
922+ assert isinstance (j11 , SubMatrixBlock )
923+
924+ assert j00 .name == 'J00'
925+ assert j01 .name == 'J01'
926+ assert j10 .name == 'J10'
927+ assert j11 .name == 'J11'
928+
929+ assert j00 .row_target == e
930+ assert j01 .row_target == e
931+ assert j10 .row_target == g
932+ assert j11 .row_target == g
933+
934+ assert j00 .col_target == e
935+ assert j01 .col_target == g
936+ assert j10 .col_target == e
937+ assert j11 .col_target == g
938+
939+ assert j00 .row_idx == 0
940+ assert j01 .row_idx == 0
941+ assert j10 .row_idx == 1
942+ assert j11 .row_idx == 1
943+
944+ assert j00 .col_idx == 0
945+ assert j01 .col_idx == 1
946+ assert j10 .col_idx == 0
947+ assert j11 .col_idx == 1
948+
949+ assert j00 .linear_idx == 0
950+ assert j01 .linear_idx == 1
951+ assert j10 .linear_idx == 2
952+ assert j11 .linear_idx == 3
953+
875954 # Check the number of submatrices
876955 assert jacobian .n_submatrices == 4
877956
0 commit comments