Pressure-gradient error?

133 views
Skip to first unread message

Nadeem Malik

unread,
Feb 10, 2023, 4:40:53 PM2/10/23
to Nek5000
Hi Neks, [changing the subject from "Integration Method"]

I have tried to implement calculation of the pressure gradient and adding it to the particle motion (ppiclfT subroutine). BUt I have found the dp/dx=0, dpdy=0. Clearly I am doing something wrong. (I assume that the pressure field is "pr" ?)

Can anyone help out here.

Here are the parts of my code that pertains to this:

[1] in *.usr:

! Pressure-gradient field
   call rzero(dpdx,n)
   call rzero(dpdy,n)
   call gradm1(dpdx,dpdy,dpdz,pr)
   call ppiclf_solve_InterpFieldUser(PPICLF_R_JDPDX,dpdx(1,1,1,1))
   call ppiclf_solve_InterpFieldUser(PPICLF_R_JDPDY,dpdy(1,1,1,1))

   if(ppiclf_npart .ne. 0) then
   do i = 1,ppiclf_npart
      dpx = ppiclf_rprop(PPICLF_R_JDPDX,i)
      dpy = ppiclf_rprop(PPICLF_R_JDPDY,i)
   end do
   end if

[2] in subroutine ppiclf_user_SetYdot:

! Pressure-gradient
fdpdx = -ppiclf_rprop(PPICLF_R_JVOLP,i)*ppiclf_rprop(PPICLF_R_JDPDX,i)
fdpdy = -ppiclf_rprop(PPICLF_R_JVOLP,i)*ppiclf_rprop(PPICLF_R_JDPDY,i)

! set ydot for all PPICLF_SLN number of equations
ppiclf_ydot(PPICLF_JVX,i) = (fbx+fqsx+fdpdx+fcx)/rmass
ppiclf_ydot(PPICLF_JVY,i) = (fby+fqsy+fdpdy+fcy)/rmass

ppiclf_ydotc(PPICLF_JVX,i) = -fqsx-fdpdx
ppiclf_ydotc(PPICLF_JVY,i) = -fqsy-fdpdy

Thanks
Nadeem   


Byeong-Cheon Kim

unread,
Feb 16, 2023, 12:08:37 PM2/16/23
to Nek5000
Hi Nadeem

Before adding the pressure gradient force with index "PPICLF_R_JDPDX", could you check  PPICLF_LRP value in  PPICLF_USER.h?
if you add the pressure gradient force term, PPICLF_LRP will be changed as 8 (in 2D) like below :

-----------------------------------------------------------------------------
#define PPICLF_LRS 4
#define PPICLF_LRP 8
#define PPICLF_LEE 1000
#define PPICLF_LEX 6
#define PPICLF_LEY 6
#define PPICLF_LRP_INT 3
#define PPICLF_LRP_PRO 3

#define PPICLF_JX 1
#define PPICLF_JY 2
#define PPICLF_JVX 3
#define PPICLF_JVY 4
#define PPICLF_R_JRHOP 1
#define PPICLF_R_JDP 2
#define PPICLF_R_JVOLP 3
#define PPICLF_R_JPHIP 4
#define PPICLF_R_JUX 5
#define PPICLF_R_JUY 6
#define PPICLF_R_JDPDX 7
#define PPICLF_R_JDPDY 8
#define PPICLF_P_JPHIP 1
#define PPICLF_P_JFX 2
#define PPICLF_P_JFY 3
-----------------------------------------------------------------------------

Second, yes "pr" is pressure field. 

I hope it will be helpful for you :)

-B.C.Kim

Nadeem Malik

unread,
Feb 16, 2023, 5:08:17 PM2/16/23
to Nek5000
Thanks B.C. It is all sorted now, as you suggested in the *.h file etc, and I am getting non-zero grad(p).

Do you know what order the accuracy of the grad(p) is the nek5000 code?

Thanks
Nadeem

Byeong-Cheon Kim

unread,
Feb 18, 2023, 6:00:07 AM2/18/23
to Nek5000
Unfortunately, I don't have any idea about the order of the grad(p)'s accuracy.

Thank you
-B.C. Kim

Philipp Schlatter

unread,
Feb 19, 2023, 1:08:00 PM2/19/23
to nek...@googlegroups.com

gradm1 is using the underlying spectral approximation to compute the derivatives, so it is high order. I assume you are using the PnPn formulation, so you don't need any mapping of the pressure.

Philipp

--
You received this message because you are subscribed to the Google Groups "Nek5000" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nek5000+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nek5000/a740d426-3865-4ffa-9967-c7b7d2fbfa91n%40googlegroups.com.
Reply all
Reply to author
Forward
0 new messages