[Nek5000-users] Calculation of Divergence Term
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Thu Nov 17 18:10:38 CST 2011
Hi,
Thanks a lot for your helpful hints and comments. I have another
problem because of boundary conditions. The special type of BC, given
in my previous message, is equivalent to no flux for the concentration
c and it must applied on all boundaries. So, the value of the integral
of the concentration over the whole volume should be conserved and
does not change with time. But after some time, c reaches very high
values which causes the Helmholtz solver to fail. In some numerical
methods, when we deal with the no flux BC, a single node on the
boundary should be excluded and be assigned a specific value. I'm not
sure if it is the case here. I wonder if you have any idea how to
implement this BC and keep the whole concentration constant.
Best Regards,
Alireza Karimi
Quoting nek5000-users at lists.mcs.anl.gov:
>
> 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
>>
> _______________________________________________
> 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
More information about the Nek5000-users
mailing list