code_saturne User's Forum

Hello,
To solve a new equations system , I want to add a new term, which depends on pressure correction (delta P), to the pressure equation (resolp)
The added term is an array called [drhodp (iel)] multiplied by the pressure correction (delta p) gives :  drhodp (iel)* (delta p).
My idea is to add directly this term to the matrix diagonal (dam) then after the "matrix" routine is called, so I can write:
do iel = 1 , ncel
  dam(iel) = dam(iel) - drhodp(iel)*volume(iel)/dt(iel)
enddo
 
Is there a better solution to consider this new term ?
 
Thanks for help.
 
Rezki
 

Hello,
Could you write your new equations system?
Best regards,
JF 

Hello,
Thanks for answering.
The system is as follows:
  * state law to resolve the density and obtain (drho/dp).
  * resolve momentum equation (preduv) using the new density.
  * resolve the new pressure equation.
Find attached the modified pressure equation.
Best regards,
RC

Rezki,
I think it is ok. The new term is fully implicit, but I am not sure that linear solvers are able to solve the new system.
Another way is to add the new term in the RHS, and add sweep in the fixed point (and update dp) to implicit it. In this case the matrix is unchanged and the system resolution is not a problem.
JF

Hi,

Thanks for your help.

Indeed, I already added this term to the RHS inside the loop, it gives :

    -------------------------------------------------------------
    do 100 isweep = 1, nswmpr
       smbr(iel) = - w7(iel) - smbr(iel)
        do iel = 1 , ncel
           smbr(iel) =  smbr(iel) + drhodp(iel)*rtp(iel,ipriph)*volume(iel)/dt(iel)
        enddo
       .....
    100  continue
----------------------------------------------------------------

Do you think I can comment the "dam"'s loop : (do iel = 1 , ncel    
dam(iel) = dam(iel) - drhodp(iel)*volume(iel)/dt(iel)       enddo) ?

Thanks in advance.

RC

hello,
You have 2 solutions:
1) add the new term in "dam" and try to solve the new system with linear solvers available in CS (it may be possible)
2) add the new term in RHS, and loop on the system to make it implicit
 
Best regards,
JF

Hello,
When trying to solve the previous system , I get divergence Jacobi error :
      "Jacobi [Pressure] : divergence après 2542 itérations : résidu initial :  8.3838e-04 ; résidu courant :  1.0069e+02"
     Pile d'appels :
         1: 0x1aba0a     <bft_error+0x3a>                 (libbft.so.1)
         2: 0xecc351     <reslin_+0x1f1>                  (libsaturne.so.0)
         3: 0x10395ed    <invers_+0x1ed>                  (libsaturne.so.0)
         4: 0x809d01c    <resple_+0x420c>                 (cs_solver)
         5: 0x80879c3    <navsto_+0x34e3>                 (cs_solver)
         6: 0x80a67c6    <tridim_+0x774a>                 (cs_solver)
         7: 0xedff02     <caltri_+0x3df2>                 (libsaturne.so.0)
         8: 0xec24f6     <cs_run+0x846>                   (libsaturne.so.0)
         9: 0xec2855     <main+0x255>                     (libsaturne.so.0)
        10: 0xb5bbd6     <__libc_start_main+0xe6>         (libc.so.6)
        11: 0x804d231    <>                               (cs_solver)
    Fin de la pile
 
The linear solver seems to be unable to inverse the pressure matrix.
Just to remind you the equation system :
1- State law , rho = f(P) using the explicit pressure, to resolve the density and obtain (drho/dp).
 
2- Resolve momentum equation (preduv) using the new density and change the "xnorme" expression :
            xnormp(iel) = xnormp(iel) +                                                                   &
                                     (propce(iel,ipcrom)-romexp(iel))/dt(iel) * volume(iel)    &
 
            romexp(iel) : Density at the previous time step.
 
3 - Resolve the new pressure equation and correct the density (attached file - see above):
 
         3-a :  w7(iel) = w7(iel) + ((propce(iel,ipcrom)-romexp(iel))/dt(iel))*volume(iel)
  
         3-b : add (drho/dp) to "dam" then "matrix" routine is called: dam(iel) = dam(iel) - drhodp(iel)*volume(iel)/dt(iel).
                  or
                  add this term to "smbr" then "itrgrp" routine is called inside isweep's loop : smbr(iel) =  smbr(iel) +  drhodp(iel)*rtp(iel,ipriph)*volume(iel)/dt(iel)
 
         3-c : Finally correct the density:
                           
  propce (iel,ipcrom) =  propce (iel,ipcrom)                                                  &
                                                            + drhodp(iel)*(rtp(iel,ipriph)-rtpa(iel,ipriph))
 
Any suggestions please ?!
 
RC

Hello,
You can try Conjugate Gradient (CG) solver or Bi-CG solver (iresol = 0 or iresol = 2).
If it still does not work, you must add the new term in the RHS.
Best regards,
JF