[Nek5000-users] Calculation of Divergence Term
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Fri Nov 18 11:18:16 CST 2011
Hi Alireza,
For the circumstances you describe, we should easily
conserve the concentration to about 8 decimal places,
so I'm guessing that there is something amiss somewhere.
If you'd like to send a tarfile with .rea/.usr and SIZE
off-list, I can take a look.
Paul
On Fri, 18 Nov 2011, nek5000-users at lists.mcs.anl.gov wrote:
> Hi,
>
> Thanks for your time. In the system which I'm simulating, the base flow is
> zero and the flow field is induced by the gradient of the concentration. So,
> Reynolds and Peclet numbers cannot be well defined. But their order of
> magnitude based on the average value of velocity is: Re~0.5 and Pe~10. So,
> the flow is not turbulent. Also, I'm using 3rd order time integration.
>
> I'm not sure about the resolving of the boundary layers. My domain is a box
> with the size of 1x1x1 and I have 2 elements in each direction (8 elements at
> all) with the degree of polynomial equal to 11. I wonder if I have to use
> finer grid with/without nonuniform spacing close to the boundaries.
>
> Regards,
> Alireza Karimi
>
>
>
>
> Quoting nek5000-users at lists.mcs.anl.gov:
>
>>
>> Hi Alireza,
>>
>> What Peclet and Reynolds numbers are you considering? How well
>> resolved are your boundary layers? Is the background flow turbulent?
>> What order of timestepping are you using (as set by TORDER)?
>> The concentration transport is not evaluated in conservation form,
>> so drift is possible, particularly if you are not fully resolved.
>>
>> I see no reason for fixing a boundary value unless you are interested in a
>> steady state solution, in which case the solution
>> is indeterminant up to a constant. For such situations in nek we
>> typically fix the mean to be zero. This should not be an issue, however,
>> for a time-marching approach to steady state.
>>
>> Let me know if you'd like me to dig deeper into the issue you are
>> having.
>>
>> Paul
>>
>>
>>
>> On Thu, 17 Nov 2011, nek5000-users at lists.mcs.anl.gov wrote:
>>
>>> 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
>>> _______________________________________________
>>> 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
> _______________________________________________
> 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