[Nek5000-users] Calculation of Divergence Term

nek5000-users at lists.mcs.anl.gov nek5000-users at lists.mcs.anl.gov
Tue Nov 15 10:30:05 CST 2011 


Alireza,

nek supports a Robin boundary condition - aka Newton law of cooling.

The standard form for this is:

   -k dT/dn = h*(T-T_inf)

where k = vdiff(.,.,.,.,ifld) is the conductivity inside the medium

      h = heat transfer coefficient

      T_inf = temperature outside the domain

As a user, you can trigger this BC by using the "c  " (convective)
boundary condition specification.  Lower case "c" indicates that
the user-supplied quantities are specified in subroutine userbc().
These quantities are:

     T_inf
     h

and passed through the NEKUSE common block.  I believe that the
variable names are "tinf" and "h" ... but am not at a terminal
where I can verify this at the moment.  The conductivity is
stored in vdiff() and may be variable or constant. 

In your case, you need something like the following.


      subroutine userbc(ix,iy,iz,iside,eg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      integer e,eg

      common /myvel/ ux(lx1,ly1,lz1,lelt),uy(lx1,ly1,lz1,lelt)
     $             , uz(lx1,ly1,lz1,lelt)

      e = gllel(eg)

      uvx = ux(ix,iy,iz,e) + vx(ix,iy,iz,e)
      uvy = uy(....
      uvz = uz ....

      if (iside.eq.1.or.iside.eq.3) then
         j1=ix
         j2=iz
      elseif (iside.eq.2.or.iside.eq.4) then
         j1=iy
         j2=iz
      else
         j1=ix
         j2=iy
      endif

      c = uvx*unx(j1,j2,iside,e)+uvy*uny(j1,j2,iside,e)+uvz*unz(j1,j2,iside,e) (this line too long)

      tinf = 0
      h = c/vdiff(ix,iy,iz,e)   ! be sure to check the sign here

      return
      end

...

hth,

Paul
      







----- Original Message -----
From: nek5000-users at lists.mcs.anl.gov
To: nek5000-users at lists.mcs.anl.gov
Sent: Monday, November 14, 2011 10:38:57 PM
Subject: Re: [Nek5000-users] Calculation of Divergence Term

Hi Alireza,

If you can cast your advection-diffusion equation for c in terms of new advection vector w = u + V then you get a standard mixed/Robin BC

dc/dn = w.n * c

that can be prescribed with existing Nek5000 infrastructure by specifying the appropriate variables from NEKUSE in userbc such as flux,hc,tinf.

Let me know if you need more details.
Aleks


----- Original Message -----
From: nek5000-users at lists.mcs.anl.gov
To: nek5000-users at lists.mcs.anl.gov
Sent: Monday, November 14, 2011 8:33:57 PM
Subject: Re: [Nek5000-users] Calculation of Divergence Term

Hi Neks,

Considering the advection-diffusion equation discussed in my previous  
messages, I want to impose a mixed type of boundary condition for the  
concentration c. In general, on the boundaries, this BC takes the  
following form as,

dc/dn = (u.n + V.n)*c

where n is the normal vector to the boundary, u is the velocity field,  
and V is the given vector field. I wonder if you have any idea how to  
impose this mixed Dirichlet-Neumann BC in the usr file.

Regards,
Alireza Karimi







Quoting nek5000-users at lists.mcs.anl.gov:

>
> No problem... sorry for the delay in getting you what you needed.
>
>
>
> On Sat, 12 Nov 2011, nek5000-users at lists.mcs.anl.gov wrote:
>
>> Hi,
>>
>> Thanks a lot for your hints and I appreciate your help. Sorry for  
>> interrupting when you were so busy.
>>
>> Regards,
>> Alireza Karimi
>>
>>
>>
>>
>> Quoting 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
>>>>>
>>>>
>>>>
>>>>
>>>> _______________________________________________
>>>> Nek5000-users mailing list
>>>> 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
>>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>>>
>>
>> _______________________________________________
>> Nek5000-users mailing list
>> 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
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>


_______________________________________________
Nek5000-users mailing list
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
https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users

More information about the Nek5000-users mailing list