[Nek5000-users] 2D Axisymmetric in semi-circular annulus
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Mon Sep 3 09:06:16 CDT 2012
Hi Aleks,
I've tested this update and it appears to have solved that problem. Thanks for looking into this!
Best wishes,
Adrian
> Hi Adrian,
>
> We have committed fixes that resolves the issues you had earlier (revision 856).
>
> Let us now if this works for you.
>
> Best.
> Aleks
>
> ----- Original Message -----
> From: nek5000-users at lists.mcs.anl.gov
> To: nek5000-users at lists.mcs.anl.gov
> Sent: Tuesday, August 28, 2012 2:19:00 AM
> Subject: Re: [Nek5000-users] 2D Axisymmetric in semi-circular annulus
>
> OK --- that does not appear to be the issue... will dig a bit
> deeper.
>
> Paul
>
> On Tue, 28 Aug 2012, nek5000-users at lists.mcs.anl.gov wrote:
>
>> Hi Adrian,
>>
>> My suspicion is that the issue is somehow related to the SYM bc + curve sides
>> + ifaxis. Generally, SYM on
>> curve sides implies using the stress formulation.
>>
>> At present, I'm not certain of the status of the combination:
>>
>> IFSTRS T
>>
>> IFAXIS T
>>
>> We can look into this as well.
>>
>> Sorry that this has been holding you up.
>>
>> Paul
>>
>> On Mon, 27 Aug 2012, nek5000-users at lists.mcs.anl.gov wrote:
>>
>>> Dear Nek users,
>>>
>>> I have been attempting to set up a 2D axisymmetric problem with an
>>> azimuthal
>>> velocity, in a semi-circular annulus container (the 2D axisymmetric version
>>> of a spherical shell), but have encountered problems with the terms
>>> evolving
>>> the azimuthal velocity blowing up near to the axis, and wondered if anybody
>>> would be able to help with this?
>>>
>>> I have looked at the vortex2 example, and have copied the relevant aspects
>>> to
>>> my .usr and .rea files. I have set IFAXIS=T to get the correct viscous
>>> operators for cylindrical geometry where (x,y)->(z,R), as well as IFHEAT=T,
>>> IFAZIV=T, using temperature as the azimuthal velocity. My test mesh is
>>> created using prenek, and contains 4 elements on a polar grid from r=0.5 to
>>> 1
>>> and theta=0 to 180. Since elements 2 & 4 have one side touching the R=0
>>> (y=0)
>>> axis, I have used "A" boundary conditions, which are applied to side 1 for
>>> each of these elements, as is required. The remaining BCs are 'SYM' and 'I'
>>> for all of the exterior sides. The .rea file containing the mesh with 4
>>> elements is attached to this email. (In the SIZE file I have set
>>> lx1=20,lxd=30,lx2=lx1-2).
>>>
>>> For a simple test, if I set the initial conditions to be a uniform rotation
>>> about the z-axis (setting temp = y), within several time-steps the code
>>> blows
>>> up, giving NaN's for each quantity (I have tested various initial
>>> conditions,
>>> which each give the same error). The line of code that appears to cause
>>> problems is the line ~240 of hmholtz.f, TERM1 =
>>> BM1(I,J,1,IEL)*U(I,J,1,IEL)/YM1(I,J,1,IEL)**2,
>>> since, when this term is set to 0, the code does not produce NaN's. (Though
>>> the viscous operator for the azimuthal velocity is then incorrect...)
>>>
>>> When IFHEAT=IFAZIV=F, i.e. when I am not integrating the
>>> temperature/azimuthal velocity, the code runs without problem on the same
>>> mesh.
>>>
>>> I am wondering if there is a problem with how I have set up the mesh that
>>> causes this problem? And if so, is there a way to fix this?
>>>
>>> Thanks,
>>> Adrian
>>>
>>> Modifications to .usr file are:
>>>
>>> c-----------------------------------------------------------------------
>>> subroutine userf (ix,iy,iz,eg)
>>> include 'SIZE'
>>> include 'TOTAL'
>>> include 'NEKUSE'
>>> integer e,f,eg
>>> if(y.gt.0) ffy = temp*temp/y
>>> return
>>> end
>>> c-----------------------------------------------------------------------
>>> subroutine userq (ix,iy,iz,eg)
>>> include 'SIZE'
>>> include 'TOTAL'
>>> include 'NEKUSE'
>>> integer e,f,eg
>>> if(y.gt.0) then
>>> visc = param(2)
>>> qvol = -uy*temp/y
>>> endif
>>> return
>>> end
>>> c-----------------------------------------------------------------------
>>> subroutine useric (ix,iy,iz,ieg)
>>> include 'SIZE'
>>> include 'TOTAL'
>>> include 'NEKUSE'
>>> temp = y
>>> return
>>> end
>>> c-----------------------------------------------------------------------
>>>
>>> My .rea file:
>>>
>>> ****** PARAMETERS *****
>>> 2.610000 NEKTON VERSION
>>> 2 DIMENSIONAL RUN
>>> 118 PARAMETERS FOLLOW
>>> 1.00000 p001 DENSITY
>>> 0.10000E-06 p002 VISCOS
>>> 0.00000 p003
>>> 1.00000 p004
>>> 0.00000 p005
>>> 0.00000 p006
>>> 1.00000 p007 RHOCP
>>> 0.100000E-06 p008 CONDUCT
>>> 0.00000 p009
>>> 0.00000 p010 FINTIME
>>> 0.100000E+08 p011 NSTEPS
>>> 0.100000E-03 p012 DT
>>> 0.00000 p013 IOCOMM
>>> 0.00000 p014 IOTIME
>>> 1000.00 p015 IOSTEP
>>> 0.00000 p016 PSSOLVER: 0=default
>>> 1.00000 p017
>>> 0.500000E-01 p018 GRID < 0 --> # cells on screen
>>> -1.00000 p019 INTYPE
>>> 5.00000 p020 NORDER
>>> 0.100000E-06 p021 DIVERGENCE
>>> 0.100000E-08 p022 HELMHOLTZ
>>> 0.00000 p023 NPSCAL
>>> 0.100000E-01 p024 TOLREL
>>> 0.100000E-01 p025 TOLABS
>>> 0.450000 p026 COURANT/NTAU
>>> 3.00000 p027 TORDER
>>> 1.00000 p028 TORDER: mesh velocity (0: p28=p27)
>>> 0.00000 p029 = magnetic visc if > 0, = -1/Rm if < 0
>>> 0.00000 p030 > 0 ==> properties set in uservp()
>>> 0.00000 p031 NPERT: #perturbation modes
>>> 0.00000 p032 #BCs in re2 file, if > 0
>>> 0.00000 p033
>>> 0.00000 p034
>>> 0.00000 p035
>>> 0.00000 p036
>>> 0.00000 p037
>>> 0.00000 p038
>>> 0.00000 p039
>>> 0.00000 p040
>>> 0.00000 p041 1-->multiplicative SEMG
>>> 0.00000 p042 0=gmres/1=pcg
>>> 0.00000 p043 0=semg/1=schwarz
>>> 0.00000 p044 0=E-based/1=A-based prec.
>>> 0.00000 p045 Relaxation factor for DTFS
>>> 0.00000 p046 reserved
>>> 0.00000 p047 vnu: mesh matieral prop.
>>> 0.00000 p048
>>> 0.00000 p049
>>> 0.00000 p050
>>> 0.00000 p051
>>> 0.00000 p052 IOHIS
>>> 0.00000 p053
>>> 0.00000 p054 fixed flow rate dir: |p54|=1,2,3=x,y,z
>>> 0.00000 p055 vol.flow rate (p54>0) or Ubar (p54<0)
>>> 0.00000 p056
>>> 0.00000 p057
>>> 0.00000 p058
>>> 0.00000 p059 !=0 --> full Jac. eval. for each el.
>>> 0.00000 p060 !=0 --> init. velocity to small nonzero
>>> 0.00000 p061
>>> 0.00000 p062 >0 --> force byte_swap for output
>>> 0.00000 p063 =8 --> force 8-byte output
>>> 0.00000 p064 =1 --> perturbation restart
>>> 0.00000 p065 #iofiles (eg, 0 or 64); <0 --> sep. dirs
>>> 0.00000 p066 output : <0=ascii, else binary
>>> 0.00000 p067 restart: <0=ascii, else binary
>>> 0.00000 p068 iastep: freq for avg_all (0=iostep)
>>> 0.00000 p069
>>> 0.00000 p070
>>> 0.00000 p071
>>> 0.00000 p072
>>> 0.00000 p073
>>> 0.00000 p074 verbose Helmholtz
>>> 0.00000 p075
>>> 0.00000 p076
>>> 0.00000 p077
>>> 0.00000 p078
>>> 0.00000 p079
>>> 0.00000 p080
>>> 0.00000 p081
>>> 0.00000 p082
>>> 0.00000 p083
>>> 0.00000 p084 !=0 --> sets initial timestep if p12>0
>>> 0.00000 p085 dt ratio if p84 !=0, for timesteps>0
>>> 0.00000 p086 reserved
>>> 0.00000 p087
>>> 0.00000 p088
>>> 0.00000 p089
>>> 0.00000 p090
>>> 0.00000 p091
>>> 0.00000 p092
>>> 20.0000 p093 Number of previous pressure solns saved
>>> 5.00000 p094 start projecting velocity after p94 step
>>> 5.00000 p095 start projecting pressure after p95 step
>>> 0.00000 p096
>>> 0.00000 p097
>>> 0.00000 p098
>>> 3.00000 p099 dealiasing: <0--> off/3--> old/4--> new
>>> 0.00000 p100
>>> 0.00000 p101 Number of additional modes to filter
>>> 1.00000 p102 Dump out divergence at each time step
>>> 0.100000E-01 p103 weight of stabilizing filter (.01)
>>> 0.00000 p104
>>> 0.00000 p105
>>> 0.00000 p106
>>> 0.00000 p107 !=0 --> add to h2 array in hlmhotz eqn
>>> 0.00000 p108
>>> 0.00000 p109
>>> 0.00000 p110
>>> 0.00000 p111
>>> 0.00000 p112
>>> 0.00000 p113
>>> 0.00000 p114
>>> 0.00000 p115
>>> 0.00000 p116 !=0: x elements for fast tensor product
>>> 0.00000 p117 !=0: y elements for fast tensor product
>>> 0.00000 p118 !=0: z elements for fast tensor product
>>> 4 Lines of passive scalar data follows2 CONDUCT; 2RHOCP
>>> 0.00000 0.00000 0.00000 0.00000 0.00000
>>> 0.00000 0.00000 0.00000 0.00000
>>> 0.00000 0.00000 0.00000 0.00000 0.00000
>>> 0.00000 0.00000 0.00000 0.00000
>>> 12 LOGICAL SWITCHES FOLLOW
>>> T IFFLOW
>>> T IFHEAT
>>> T IFTRAN
>>> T T F F F F F F F F F IFNAV & IFADVC (convection in P.S. fields)
>>> F F T T T T T T T T T T IFTMSH (IF mesh for this field is T mesh)
>>> T IFAXIS
>>> T IFAZIV
>>> T IFSTRS
>>> F IFMGRID
>>> F IFMODEL
>>> F IFKEPS
>>> F IFCHAR
>>> 10.0000 10.0000 -4.00000 -5.50000
>>> XFAC,YFAC,XZERO,YZERO
>>> **MESH DATA** 1st line is X of corner 1,2,3,4. 2nd line is Y.
>>> 4 2 4 NEL,NDIM,NELV
>>> ELEMENT 1 [ 1C] GROUP 0
>>> -0.4371138E-07 -0.2185569E-07 0.3535533 0.7071066
>>> 1.000000 0.5000000 0.3535533 0.7071066
>>> ELEMENT 2 [ 1D] GROUP 0
>>> 0.5000000 1.000000 0.7071066 0.3535533
>>> 0.000000 0.000000 0.7071066 0.3535533
>>> ELEMENT 3 [ 1D] GROUP 0
>>> 0.5962440E-08 0.1192488E-07 -0.7071066 -0.3535533
>>> 0.5000000 1.000000 0.7071066 0.3535533
>>> ELEMENT 4 [ 1 ] GROUP 0
>>> -1.000000 -0.5000000 -0.3535533 -0.7071066
>>> 0.8742278E-07 0.4371139E-07 0.3535533 0.7071066
>>> ***** CURVED SIDE DATA *****
>>> 8 Curved sides follow IEDGE,IEL,CURVE(I),I=1,5, CCURVE
>>> 2 1 -0.500000 0.00000 0.00000 0.00000 0.00000
>>> C
>>> 4 1 1.00000 0.00000 0.00000 0.00000 0.00000
>>> C
>>> 2 2 1.00000 0.00000 0.00000 0.00000 0.00000
>>> C
>>> 4 2 -0.500000 0.00000 0.00000 0.00000 0.00000
>>> C
>>> 2 3 1.00000 0.00000 0.00000 0.00000 0.00000
>>> C
>>> 4 3 -0.500000 0.00000 0.00000 0.00000 0.00000
>>> C
>>> 2 4 -0.500000 0.00000 0.00000 0.00000 0.00000
>>> C
>>> 4 4 1.00000 0.00000 0.00000 0.00000 0.00000
>>> C
>>> ***** BOUNDARY CONDITIONS *****
>>> ***** FLUID BOUNDARY CONDITIONS *****
>>> E 1 1 3.00000 1.00000 0.00000 0.00000 0.00000
>>> SYM 1 2 0.00000 0.00000 0.00000 0.00000 0.00000
>>> E 1 3 2.00000 3.00000 0.00000 0.00000 0.00000
>>> SYM 1 4 0.00000 0.00000 0.00000 0.00000 0.00000
>>> A 2 1 0.00000 0.00000 0.00000 0.00000 0.00000
>>> SYM 2 2 0.00000 0.00000 0.00000 0.00000 0.00000
>>> E 2 3 1.00000 3.00000 0.00000 0.00000 0.00000
>>> SYM 2 4 0.00000 0.00000 0.00000 0.00000 0.00000
>>> E 3 1 1.00000 1.00000 0.00000 0.00000 0.00000
>>> SYM 3 2 0.00000 0.00000 0.00000 0.00000 0.00000
>>> E 3 3 4.00000 3.00000 0.00000 0.00000 0.00000
>>> SYM 3 4 0.00000 0.00000 0.00000 0.00000 0.00000
>>> A 4 1 0.00000 0.00000 0.00000 0.00000 0.00000
>>> SYM 4 2 0.00000 0.00000 0.00000 0.00000 0.00000
>>> E 4 3 3.00000 3.00000 0.00000 0.00000 0.00000
>>> SYM 4 4 0.00000 0.00000 0.00000 0.00000 0.00000
>>> ***** THERMAL BOUNDARY CONDITIONS *****
>>> E 1 1 3.00000 1.00000 0.00000 0.00000 0.00000
>>> I 1 2 0.00000 0.00000 0.00000 0.00000 0.00000
>>> E 1 3 2.00000 3.00000 0.00000 0.00000 0.00000
>>> I 1 4 0.00000 0.00000 0.00000 0.00000 0.00000
>>> A 2 1 0.00000 0.00000 0.00000 0.00000 0.00000
>>> I 2 2 0.00000 0.00000 0.00000 0.00000 0.00000
>>> E 2 3 1.00000 3.00000 0.00000 0.00000 0.00000
>>> I 2 4 0.00000 0.00000 0.00000 0.00000 0.00000
>>> E 3 1 1.00000 1.00000 0.00000 0.00000 0.00000
>>> I 3 2 0.00000 0.00000 0.00000 0.00000 0.00000
>>> E 3 3 4.00000 3.00000 0.00000 0.00000 0.00000
>>> I 3 4 0.00000 0.00000 0.00000 0.00000 0.00000
>>> A 4 1 0.00000 0.00000 0.00000 0.00000 0.00000
>>> I 4 2 0.00000 0.00000 0.00000 0.00000 0.00000
>>> E 4 3 3.00000 3.00000 0.00000 0.00000 0.00000
>>> I 4 4 0.00000 0.00000 0.00000 0.00000 0.00000
>>> 0 PRESOLVE/RESTART OPTIONS *****
>>> 7 INITIAL CONDITIONS *****
>>> C Default
>>> C Default
>>> C Default
>>> C Default
>>> C Default
>>> C Default
>>> C Default
>>> ***** DRIVE FORCE DATA ***** BODY FORCE, FLOW, Q
>>> 4 Lines of Drive force data follow
>>> C
>>> C
>>> C
>>> C
>>> ***** Variable Property Data ***** Overrrides Parameter data.
>>> 1 Lines follow.
>>> 0 PACKETS OF DATA FOLLOW
>>> ***** HISTORY AND INTEGRAL DATA *****
>>> 0 POINTS. Hcode, I,J,H,IEL
>>> ***** OUTPUT FIELD SPECIFICATION *****
>>> 6 SPECIFICATIONS FOLLOW
>>> T COORDINATES
>>> T VELOCITY
>>> T PRESSURE
>>> T TEMPERATURE
>>> F TEMPERATURE GRADIENT
>>> 0 PASSIVE SCALARS
>>> ***** OBJECT SPECIFICATION *****
>>> 0 Surface Objects
>>> 0 Volume Objects
>>> 0 Edge Objects
>>> 0 Point Objects
More information about the Nek5000-users
mailing list