[Nek5000-users] Calculation of Divergence Term

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


Yes... sorry... you have to change the names.

Thanks,

Paul

PS - also, noted Alek's earlier email - I think "hc" is the ht trnsfer 
coef.



On Tue, 15 Nov 2011, nek5000-users at lists.mcs.anl.gov wrote:

> Hi Paul,
>
> If I'm not mistaken, aren't the variable names ux, uy, and uz already
> taken in the common blocks in NEKUSE? If so, wouldn't the declaration
> of the myvel common block violate those names? Just wanted to make
> sure Alireza does not run into any pesky compilation errors that may
> be difficult to spot.
>
> Josh
>
>
>
> On Nov 15, 2011, at 9:30 AM, "nek5000-users at lists.mcs.anl.gov"
> <nek5000-users at lists.mcs.anl.gov> wrote:
>
>>
>>
>> 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
>> _______________________________________________
>> 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