[Nek5000-users] Questions in AXHELM and WLAPLACIAN

Sat May 31 18:19:50 CDT 2014

Hi Emmanuel,

Regarding the gradm1() crash, did you declare Tx,Ty, and Tz in a new common block
of your own creation?  This would be important to do.

Regarding where to modify the code - the most common approach is to create a new
routine (say, "subroutine make_term3") which generates the required source terms
when called once per timestep from userchk.   You would store these source terms
in an array and then load from this array in the userq() routine, which is called
point-by-point, field-by-field, by nek as needed.

Regarding your earlier email about replacing axhelm, the strategy you
propose would not give the desired result.   My suggestion would be
to look carefully at some of the FEM/SEM literature.   There are a couple
of JCP papers by Tomboulides & Orszag (for low Mach), my 97 JCP paper
(on Schwarz, but generally on the SEM discretization), and the book
High Order Methods for Incompressible Flow with Deville & Mund.
Please let me know if you have further questions.

Best, Paul

________________________________
From: nek5000-users-bounces at lists.mcs.anl.gov [nek5000-users-bounces at lists.mcs.anl.gov] on behalf of nek5000-users at lists.mcs.anl.gov [nek5000-users at lists.mcs.anl.gov]
Sent: Saturday, May 31, 2014 3:47 AM
To: nek5000-users at lists.mcs.anl.gov
Subject: Re: [Nek5000-users] Questions in AXHELM and WLAPLACIAN

Hi Paul

I have some difficulties to compute a simple gradient of a variable. I think I’m confusing something.

Basically I would like to modify the temperature and species equations to add a new term in order to take into account a velocity correction related to mass conservation and species diffusion coefficients.

For temperature the equation is:

d T /dt = term1 + term2 + term3
with
term1 = - U grad (T)
term2 =   1/rho Cp  * grad ( lambda * grad (T))
term3 =   1/rho Cp  * (rho * sum_k [Cp,k * (Y_k * D_k * grad(Y_k) - Y_k * sum_j [D_j * grad(Y_j )]  ) ] ) * grad(T)

and for species equation :

d Y /dt = term1 + term2 + term3
with
term1 = - U grad (Y)
term2 =   1/rho Cp  * grad ( lambda * grad (Y))
term3 = - 1/rho * grad (Y * sum_k [D_k * grad (Y_k) ] )

where D_k or D_j are the species diffusion coefficients, T the temperature and Y the species mass fraction.

I am using the low-Mach number formulation and solving the equations with CVODE.
Basically term1 is computed by the routine convab while term2 is computed by routine wlapacian. This is the current official nek implementation.

Now I would like to add the term3, which involves some gradients of temperature or species variables. I think the simplest way is to modify the routine convab or maybe makeq.f

However I tried a simple call gradm1 (Tx,Ty,Tz,T) but NEK crashes with a segmentation fault error. I tried to debug the routine gradm1, and it appears that it crashes during the loop

do i=1,lxyz
               ux(i,e) =jacmi(i,e)*(ur(i)*rxm1(i,1,1,e)
     $                            + us(i)*sxm1(i,1,1,e) )
               uy(i,e) =jacmi(i,e)*(ur(i)*rym1(i,1,1,e)
     $                            + us(i)*sym1(i,1,1,e) )
            enddo

after reaching a particular element e (in my case it crashes for element 9).

So my questions are just:
1) According to the implementation of term3 in temperature and species equations, what is the best way to proceed? I have think about modifying convab or makeq but is it a good idea?

2) Is it a good idea to use gradm1? Why does it crash for a particular element? Should I pay attention to something in particular when using this routine?

Thank you very much for your help !

Best regards,
Emmanuel

On 30 May 2014, at 10:25, nek5000-users at lists.mcs.anl.gov<mailto:nek5000-users at lists.mcs.anl.gov> wrote:

Hi Emmanuel,

Strictly speaking, our functions are only once differentiable -- they are C0, not C1.

The test functions are also in H1 (i.e., are C0) and therefor the standard Galerkin
approach is to transfer one derivative to the test functions via integration by parts
and thus avoid twice differentiating a function.

If you just want the gradient of something, you would use gradm1 (gradient for
a quantity on mesh 1).

Paul

________________________________
From: nek5000-users-bounces at lists.mcs.anl.gov<mailto:nek5000-users-bounces at lists.mcs.anl.gov> [nek5000-users-bounces at lists.mcs.anl.gov<mailto:nek5000-users-bounces at lists.mcs.anl.gov>] on behalf of nek5000-users at lists.mcs.anl.gov<mailto:nek5000-users at lists.mcs.anl.gov> [nek5000-users at lists.mcs.anl.gov<mailto:nek5000-users at lists.mcs.anl.gov>]
Sent: Thursday, May 29, 2014 7:47 PM
To: nek5000-users at lists.mcs.anl.gov<mailto:nek5000-users at lists.mcs.anl.gov>
Subject: Re: [Nek5000-users] Questions in AXHELM and WLAPLACIAN

Hi Paul,

Thank you for your quick answer, this is really helpful.

I still have a question concerning the computation of a gradient.

Why did you construct the diffusion operator with the helmholtz routine? Why not using directly the routines gradm1 ?

By the way I am not sure to understand the differences between gradm1 (navier5.f), wgradm1 (navier1.f) , gradm11 (navier2.f) and gradm11ts (navier2.f).
I guess that the w letter refers to the weak formulation, but I don’t understand in which cases we need to choose a formulation instead of others.

In summary, if I want to modify the temperature/passive scalar equations to add a term that contains a gradient of the solved temperature/PS variable, for example, which “gradm” routine do I need to use?

Thank you for your help !

Best regards,
Emmanuel

On 29 May 2014, at 20:20, nek5000-users at lists.mcs.anl.gov<mailto:nek5000-users at lists.mcs.anl.gov> wrote:

Hi Emmanuel,

axhelm effects the weak Laplacian plus the mass matrix.

Specifically, if you wanted to solve:

(1)   -nabla . [ h1 * nabla ] u + h2 * u   =  f

with h1, h2, and f all functions of X=(x,y,z),
the discrete form would be:

(2)         A(h1) u + B(h2) u = B f

or

(3)         H u = B f

where A is the (h1-dependent stiffness matrix and
B(h2) is the h2-dependent mass matrix.

For the SEM, B is diagonal, so B(h2) = h2*B, where
we view h2 as a diagonal matrix.

In Nek5000, a matrix-vector product  w = Hu is effected via

     call axhelm (w,u,h1,h2,imesh,isd)

with w,u,h1,h2 being arrarys.   Here, imesh indicates
if you are on the fluid mesh or thermal mesh.  The fluid
mesh is a subset of the thermal mesh in the case of conjugate heat transfer between a fluid and a solid.
isd indicates which spatial direction is considered
in axisymmetric applications where u is a component
of a vector field.

Notice the sign of the Laplacian in Eq. (1) -- that may
help resolve your second question?  (There is a sign
change that ariss from the standard integration-by-parts
stepin the weak form.)

Paul

On Thu, 29 May 2014, nek5000-users at lists.mcs.anl.gov<mailto:nek5000-users at lists.mcs.anl.gov> wrote:

Dear nek5000 team,

Could you help me with the 2 questions below? Thank you very much !

1) Concerning the routine “axhelm” in hmholtz.f. If I am right, one
possibility of this routine can be the discretisation of the diffusion
operator. However it is not clear to me if axhelm performs:

a) vdiff  x   Laplacian(phi)

or

b) grad ( vdiff x grad (phi) )

phi being the variable solved for (T or PS for example).

So, is it actually a) or b) ?

2) In the routine wlaplacian in navier1.f, I am not sure to understand the
consequence of the line:

call sub2 (out,wrk,ntot)

after the lines:

call bcneusc(out,1) call axhelm(wrk,a,diff,h2,imesh,1)

If I am right, this will change the sign in front of the diffusion operator,
right?  Because in makeq.f, I am not sure that the line call
add2(bq(1,1,1,1,ifield-1),w1,ntot)

will actually assemble the RHS of the transport equation :  - (vtrans) x  U
grad (phi)  +  grad ( vdiff x grad (phi) )

Thank you for your help !

Best regards, Emmanuel

_______________________________________________
Nek5000-users mailing list
Nek5000-users at lists.mcs.anl.gov<mailto:Nek5000-users at lists.mcs.anl.gov>
https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users

_______________________________________________
Nek5000-users mailing list
Nek5000-users at lists.mcs.anl.gov<mailto:Nek5000-users at lists.mcs.anl.gov>
https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20140531/ae281fca/attachment-0001.html>