Integral COnstraint for single Solver

Numerical methods and mathematical models of Elmer
Juha
Site Admin
Posts: 357
Joined: 21 Aug 2009, 15:11

Re: Integral COnstraint for single Solver

Post by Juha »

Hi,

there is also a second system if you have constraints such that you can attach the new rows
to DOFs of the original matrix, namely A % ConstraintMatrix. An example of using
this to constraint values on a boundary to given function would be:

p => Solver % Variable % Perm
CM => AllocateMatrix()
CM % Format = MATRIX_LIST
ALLOCATE(temp(A % NumberOfRows), Indexes(Solver % MaxElementDOFs)); temp=0
DO t=1,Mesh % NumberOfBoundaryElements
Element=>ActiveBoundaryElement(t)
BC=>GetBC()
val(1:n) = GetReal(BC,'some setting',Found)
IF(.NOT.Found) CYCLE
n = GetElementNOFNodes()
nd = GetElementDOFs(Indexes)
CALL GetMassMatrix(Element,n,nd,val,MASS,FORCE) ! the usual elementwise
! integration, your own routine
DO i=1,nd
DO j=1,nd
CALL AddtoMatrixElement(CM,p(Indexes(i)),p(Indexes(j),MASS(i,j))
END DO
temp(p(Indexes(i)))=temp(p(Indexes(i))) + FORCE(i)
END DO
END DO
CALL List_toCRSMatrix(CM)
ALLOCATE(CM % Rhs(CM % NumberOfRows)); CM % Rhs=Temp(1:CM % NumberOFRows)
ALLOCATE(CM % InvPerm(CM % NumberOfRows))
DO i=1,CM % NumberofRows
CM % InvPerm(i) = i
END DO
A % ConstraintMatrix => CM

Note that you don't need to set the transpose, this is done automatically. Note also, that the CM matrix would most probably have a large number of zero rows, and one could (should?) compact the matrix using the permutation "InvPerm", not done here though.

The value of this code vs the addmatrix stuff is that the global matrix rows are only
shared by partitions as needed, and could be more scalable.

The code could be simplified if one could use DefaultUpdateEqations and friends
to update the CM matrix, unfortunately this is not really possible at the moment.
Something for the future perhaps.

The code is also completely unstested....

Regards,Juha
Franz Pichler
Posts: 196
Joined: 29 Sep 2011, 12:25
Antispam: Yes

Re: Integral COnstraint for single Solver

Post by Franz Pichler »

HI there,

i have started working on the integral constraint again and i was hoping that someone could dive into the INvPerm for me a little bit. As i understand, in the pseudo code above there are ase many rows as in the original matrix. I would like to create a constraint matrix with just one row now. So i only add to the (1,perm(indexes(j)) entry and what do i do with the invperm then?

best regards

Franz
Franz Pichler
Posts: 196
Joined: 29 Sep 2011, 12:25
Antispam: Yes

Re: Integral COnstraint for single Solver

Post by Franz Pichler »

Okay also without using the invperm i managed to make this work sor a serial run. In parallel it doesnt do the same thing though. It runs through but somehow the integral condition has a wrong effect. so i add the entries for all partitions in row number one and i add to rhs entry number one. i am realizing an integral condition so the processors should add up the entries in the right manner. is that also true for the rhs?

best regards

Franz
Post Reply