Integral COnstraint for single Solver

Numerical methods and mathematical models of Elmer
Franz Pichler
Posts: 196
Joined: 29 Sep 2011, 12:25
Antispam: Yes

Integral COnstraint for single Solver

Post by Franz Pichler »

Hi there i don't remember quite exactly, is it possible to switch on some integral conditions for a solver?
Best would be one where i can program the matrix structure by myself.
i don't find much about that in the forum (searched for constrained,restriction and integral) .

I found in the Solver manual that there are the constraints for a coupled solver. mayybe i should go for that. but i rather would use my own solver because i do alot of tweaking in the nonlinear loop. and as i understand the coupled solvers use pretty much only the bulkassembly routines of the solver. is that right?

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 »

Hi there,
just a little addon:

i don't even need a real integral constraint. if i could just add a row to the system matrix where i can realize a constraint i would be happy. i guess one way to do this would be the "before linear solve" feature right?
to do irt in the solver is probaly not that easy?

best regards
Franz
raback
Site Admin
Posts: 4812
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

Re: Integral COnstraint for single Solver

Post by raback »

Hi

I think you could use some ListMatrix operations if you're not critical on computational cost. Assume you have a standard CRS matrix. You could use pseudocode like this:

Code: Select all

A=>Solver % Matrix
CALL List_ToListMatrix( A )
DO your stuff
  CALL AddToMatrixElement( A, i, j, val )
END DO
CALL List_ToCRSMatrix( A )
The list matrix adds more flexibility since the (i,j) entries do no need to match the existing entries. Once you have done this once you don't need to redo it as the matrix will have all the needed entries to set the constraints.

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

Re: Integral COnstraint for single Solver

Post by Franz Pichler »

thanks peter,

thats perfect.

could you please also give me some pseudo code to change the rhs in the same manner?

thanks in advance
have a good one ,
Franz
Franz Pichler
Posts: 196
Joined: 29 Sep 2011, 12:25
Antispam: Yes

Re: Integral COnstraint for single Solver

Post by Franz Pichler »

oh and also, does this work in parallel?
raback
Site Admin
Posts: 4812
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

Re: Integral COnstraint for single Solver

Post by raback »

Hi Franz

To my understanding, it does not work in parallel except if all the new connections are in one partition. So if you want to do some stuff for boundaries you can use "-connect bc_id" to enforce that all elements in given boundary are in the same partition. This works best with Metis & dual graph e.g. "-metis Np 3 -partdual".

I would do the loop once as a dummy before going for the actual physical assembly. However, I do see your point in that the r.h.s. values need to be enlarged as well. Of course this depends on the nature of your constraint. Some constraints don't add additional lines of code.

There might be a better strategy if you need to add entries. Just create an additional matrix "Matrix % CM" (as in constraint matrix) and this will be incorporated into the system when solving the primary system. Unfortunately there aren't many user examples but somewhere in the code there could be pieces. Below is a mock-up that I simplified from some recent code.

Code: Select all

  PS => Asolver % Variable % Perm
  nm =  Asolver % Matrix % NumberOfRows
  IF(.NOT.ASSOCIATED(CM)) THEN
    CM => AllocateMatrix()
    CM % FORMAT = MATRIX_LIST
    ALLOCATE(CM % RHS(nm + 1)); CM % RHS=0._dp
  ELSE
    CM % RHS=0._dp
    IF(ASSOCIATED(CM % Values)) CM % Values=0._dp;
  END IF

  Prm => Solver % Variable % Perm
  IF (ParEnv % myPE == 0) THEN
    lambda = ...
    CM % RHS(nm+1) = lambda 
    CALL SetMatrixElement(CM,nm+1,nm+1,0.0)
  END IF

  ! Rather loop just over active ones, this is just demo  
   DO i=1,nm
      coeff = ....
      CALL AddToMatrixElement(CM,nm+1,i,coeff)
   END DO

  IF(CM % FORMAT==MATRIX_LIST) CALL List_toCRSMatrix(CM)
  Asolver %  Matrix % AddMatrix => CM
-Peter

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

Re: Integral COnstraint for single Solver

Post by Franz Pichler »

Hi Peter thanks for this hint,

that looks promising.

Do i need the global indexes of the nodes in the AddToMatrixElement here?

so if i wanted to constraint an element such that the average value of a variable on the elements nodes is say constraint_value.
can i just add 1/numberof elementnodes at Element%Nodeindex(i) for every i?
or do i have to find out the global index in the matrix? ( i am thinkin of parallel again)

bye
Franz
raback
Site Admin
Posts: 4812
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

Re: Integral COnstraint for single Solver

Post by raback »

Hi

The global numbering should be taken care of by the system.

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

Re: Integral COnstraint for single Solver

Post by Franz Pichler »

HI Peter:
i am still struggling a bit.

my curretn flow is the following,

in the assembly loop i have

Code: Select all

                     CALL DefaultInitialize()
                     DO t = 1, Solver % NumberOfActiveElements
                              Element => GetActiveElement(t)
                              n  = GetElementNOFNodes(Element)
                              nd = GetElementNOFDOFs(Element)
                              IF( GetElementFamily() == 1 ) CYCLE
                              CALL LocalMatrix( Element, n, nd)
                      end do  
                      print *,"ASSEMBLED"

                       Element => GetActiveElement(ElementToConstraint)
                      CALL AddConstraints( Element,constraint_value)
                      CALL DefaultFinishBulkAssembly()
                      CALL DefaultFinishAssembly()                      
                      CALL DefaultDirichletBCs()   
                      Norm = DefaultSolve()
ANd the addconstraint function is

Code: Select all


  SUBROUTINE AddConstraints( Element,constraint_value)
        TYPE(Element_t) :: Element
        TYPE(Matrix_T),POINTER::CM
        REAL (KIND=dp)::constraint_value,coeff
        INTEGER :: nm
        SAVE CM

       nm =  Solver % Matrix % NumberOfRows
       IF(.NOT.ASSOCIATED(CM)) THEN
         CM => AllocateMatrix()
         CM % FORMAT = MATRIX_LIST
         ALLOCATE(CM % RHS(nm + number_of_constraints)); CM % RHS=0._dp
       ELSE 
         CM % RHS=0._dp
         IF(ASSOCIATED(CM % Values)) CM % Values=0._dp
       END IF


         CM % RHS(nm+index_of_constraint) = constraint_value
         print *,parenv%mype,nm,MODEL%NumberofNodes,number_of_constraints,constraint_value
         CALL SetMatrixElement(CM,nm+index_of_constraint,nm+constraint_i,0.0)
     
        DO i=1,SIZE(Element % NodeIndexes)
           coeff = 1./SIZE(Element % NodeIndexes)
           print *,"COEFF:",parenv%mype,i,coeff,Element % NodeIndexes(i),Solver % Variable % Perm(Element % NodeIndexes(i))
           CALL AddToMatrixElement(CM,nm+1,Solver % Variable % Perm(Element % NodeIndexes(i)),coeff)
        END DO
      
       IF(CM % FORMAT==MATRIX_LIST) CALL List_toCRSMatrix(CM)
       Solver %  Matrix % AddMatrix => CM
    END SUBROUTINE AddCOnstraints

SO my first question is:
How are the "defaultfinish..." and the "addconstraint " to arrange in the upper part?

ANd i did not quite understand the check if parenv%mype ==0 in the code that you posted above.
as i understand every proc has its own addmatrix. so if i let every proc add one constraint i have as many additional lines in the global system matrix as i have processors, right?
and so every processor should take care of its rhs?
or are the addmatrices overlapped somehow in the global(with global i mean added up of the parallel processes) system matrix?

a full example code would be great,
anyway as always all hints are highly appreciated!

an elmer fan
Franz
Juha
Site Admin
Posts: 357
Joined: 21 Aug 2009, 15:11

Re: Integral COnstraint for single Solver

Post by Juha »

Hi,
The AddMatrix system adds new global row(s) to the matrix that are shared by all partitions. You can think of them having globally indices from 1 to max_partitions(addmatrix%numberofrows-a%numberofrows). So, no, the partition AddMatrices are not independent, but share the row numering. If you need separate constraints from different partitions you need to setup the numbering by yourself (and add dummy (diagonal) entries to AddMatrix for those partitions that don't use rows numbered
lower to those they really do use, e.g. CALL AddToMatrixElement(CM,row,row,0._dp))

Note that locally you need to assign the new rows to AddMatrix starting from A % numbeofrows+1, A being the original local system matrix (this allows to modify the
original system matrix also). For Lagrange coefficient constraint system you also
add the correspoding columns (e.g. add transpose of your new equations to the
matrix). If you don't (the equations are not directly Lagrange coefficient constrains,
but some other stuff), please add zeros to the Diagonal of AddMatrix(1..A % NumberORows)

A pseudo code constraining the sum of all dofs in A to 0 would be

Parallel = ParEnv % PEs>1
me = ParEnv % myPE

CM => AllocateMatrix()
CM % Format = MATRIX_LIST
n = A % NumberOfRows
ALLOCATE(CM % Rhs(1)); CM % Rhs=0._dp
DO i=1,n
IF(Parallel) THEN
IF(A % ParallelInfo % NeighbourList(i)%Neighbours(1)/=me)CYCLE
END IF
CALL AddToMatrixElement( CM, i, n+1, 1._dp)
CALL AddToMatrixElement( CM, n+1,i, 1._dp)
END DO
IF(CM % Format==MATRIX_LIST) CALL List_toCRSMatrix(CM)
A % AddMatrix => CM

So this code adds one row & column with all entries equal to 1 (but the diagonal) to the system matrix. This should constraint the sum of all values to 0 (or the value given to CM % Rhs(1)).

The code should work both in serial and in parallel. The complication with the parallel
case is because the A matrix has duplicate rows/columns for the partition interfaces
and this is sum, not integral (the code checks that it updates the AddMatrix only from one of the partitions sharing the DOFs).

For integral constraints you should'nt have to worry 'bout sharings as the values really should be summed together from all partions sharing the DOFs.

I'd recommed updating your ElmerSolver from SVN, although this should work out, if
you didn't update yesterday or this morning ;-)

Please ask, if i missed something...
Regards, Juha
Post Reply