[Nek5000-users] Calculation of Divergence Term

[nek5000-users at lists.mcs.anl.gov]

Dear Alireza,

You would add this to qvol, since it is the convection-diffusion
equation -- so, the routine is userq.

You should precompute the term in userchk, store the result in
a common block, and then dereference in userq.   Expanding on 
Josh's suggestion below you would have something like

c-----------------------------------------------------------------------
       subroutine userchk
       :
       :

       common /mydiv/ divcv(lx1,ly1,lz1,lelt)
       common /scrns/ dcv1dx(lx1,ly1,lz1,lelt)
                    , dcv2dy(lx1,ly1,lz1,lelt)
                    , dcv3dz(lx1,ly1,lz1,lelt)
                    , work1 (lx1,ly1,lz1,lelt)
                    , work2 (lx1,ly1,lz1,lelt)


         then his calls... be sure to define n

       return
       end
c-----------------------------------------------------------------------
       subroutine userq  (ix,iy,iz,eg)
       include 'SIZE'
       include 'TOTAL'
       include 'NEKUSE'

       common /mydiv/ divcv(lx1,ly1,lz1,lelt)

       integer e,f,eg
       e = gllel(eg)

       qvol = divcv(ix,iy,iz,e)  Be sure to check the sign of this term

       return
       end
c-----------------------------------------------------------------------



It's not a priori clear that you need to dealias this term.
Dealiasing is needed in certain circumstance for stability -- it
does not materially affect accuracy.  Without further context
I cannot tell you whether you need to dealias or not.

The approach suggested by Josh is perfectly adequate for the
standard (non-dealiased) approach.


Paul


On Fri, 11 Nov 2011, nek5000-users at lists.mcs.anl.gov wrote:

> Hi,
>
> I'm sending this message again, because I got no answer for my last email. I 
> have some questions regarding the additional term given in my previous 
> message. I don't know where and in which subroutine I should add this term to 
> the right hand side of the advection-diffusion equation. Also, could you 
> please provide me with the best option I can take to dealias this term (if 
> needed at all)? I wonder if you could help me in this regard.
>
> Best Regards,
> Alireza Karimi
>
>
>
> Quoting nek5000-users at lists.mcs.anl.gov:
>
>> Alireza,
>> 
>> For this, I'll assume that you have split V into V1, V2, V3, and that
>> they have the structure
>> 
>> V1(lx1,ly1,lz1,lelv)  (and likewise for the other two)
>> 
>> 
>> First, to multiply c*V1 (V2, V3), you can use col2.  It would look like
>>
>>  call col2(V1,c,n)
>>  call col2(V2,c,n)
>>  call col2(V3,c,n)
>> 
>> where n = nx1*ny1*nz1*nelv.  A couple of things to note here:
>>
>>  1) col2 will overwrite V1 (V2, V3).  Effectively, the operation is
>> V1 = c*V1.  If this is not desirable, then you could use col3 with
>> something like
>>
>>         call col3(cV1,V1,c,n)
>>
>>     This will store it in cV1 and leave V1 untouched.
>>
>>  2) You could run into aliasing issues by multiplying c and V
>> directly, and so it may be wise to dealias the product before taking
>> the derivative.  Developers, maybe you could suggest the best way of
>> going about doing this in Nek???
>> 
>> Moving on, to take the derivatives, there are a few choices.  One way
>> would be do to something like
>>
>>   call gradm1(dcV1dx,work1  ,work2  ,cV1)
>>   call gradm1(work1  ,dcV2dy,work2  ,cV2)
>>   call gradm1(work1  ,work2  ,dcV3dz,cV3)
>>
>>   call add4(divcV, dcV1dx, dcV2dy, dcV3dz, n)
>> 
>> which would give you div(cV) in the array divcV.  Note that dcV1dx,
>> dcV2dy, dcV3dz, work1, work2, and divcV all should have the same
>> storage structure as V1, c, etc.
>> 
>> This is probably the simplest way that I can think to do it, although
>> it is somewhat inefficient since we are computing, for instance,
>> d(cV1)/dy unnecessarily.  To make it more efficient, you can manually
>> compute the derivatives individually (take a look at the gradm1
>> subroutine to get an idea of how to do this, located in navier5.f in
>> the source)
>> 
>> Hope this helps, and developers, please correct me if I did something
>> horribly wrong  here :).
>> 
>> Josh
>> 
>> On Fri, Nov 4, 2011 at 4:59 PM,  <nek5000-users at lists.mcs.anl.gov> wrote:
>>> Hi,
>>> 
>>> I want to solve an advection-diffusion equation in conjunction with the
>>> momentum and mass conservation equations. But my equation has an extra
>>> divergence term which I don't know how to calculate it in Nekton. So,
>>> denoting the scalar concentration field with c and velocity field with u,
>>> the equation turns out to be:
>>> 
>>> dc/dt + u.grad c = Laplacian c - div (c V)
>>> 
>>> which by d/dt I mean partial derivative with respect to time. V is a known
>>> vector field calculated from the vorticity tensor. I wonder if you could
>>> show me how to evaluate this additional term (the last term) in the user
>>> file.
>>> 
>>> 
>>> Best Regards,
>>> Alireza Karimi
>>> 
>>> 
>>> _______________________________________________
>>> Nek5000-users mailing list
>>> Nek5000-users at lists.mcs.anl.gov
>>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>>> 
>> 
>> 
>> 
>> --
>> Josh Camp
>> 
>> "All that is necessary for the triumph of evil is that good men do
>> nothing" -- Edmund Burke
>> _______________________________________________
>> Nek5000-users mailing list
>> Nek5000-users at lists.mcs.anl.gov
>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>> 
>
>
>
> -- 
> Alireza Karimi
> PhD Candidate
> Department of Engineering Science and Mechanics
> Virginia Polytechnic Institute and State University
> Blacksburg, VA 24061
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users