[Nek5000-users] 2D Axisymmetric in semi-circular annulus
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Fri Aug 31 22:08:56 CDT 2012
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
>> _______________________________________________
>> 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