Howdy Neks, is there a function in Nek that for a given point in x,y,z would return processor number, element number and r,s,t index of the nearest GLL grid point? I was desperately trying to figure out how the interpolation routines (f.e. for history points) would help me, but I only got garbage (-1, 0 and 1) back for the variables I would expect contain r,s,t. It worked with the element number, however. Thanks, Markus
Hi Markus, indeed hpts() is a good example to look at. The array rst contains the canonical coordinates you're looking for. Note the rst is a packed array (r_p1, s_p1, t_p1,r_p2,s_p2,t_p2, ...). Cheers, Stefan On Jul 9, 2010, at 9:46 AM, [email protected] wrote:
Howdy Neks,
is there a function in Nek that for a given point in x,y,z would return processor number, element number and r,s,t index of the nearest GLL grid point? I was desperately trying to figure out how the interpolation routines (f.e. for history points) would help me, but I only got garbage (-1, 0 and 1) back for the variables I would expect contain r,s,t. It worked with the element number, however.
Thanks, Markus
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Markus, those r,s,t values (not indices) are valid - do they not correspond to your expectations? Off-hand one might expect that the associated indices would be 1,NY1/2,NZ1 Paul On Fri, 9 Jul 2010, [email protected] wrote:
Howdy Neks,
is there a function in Nek that for a given point in x,y,z would return processor number, element number and r,s,t index of the nearest GLL grid point? I was desperately trying to figure out how the interpolation routines (f.e. for history points) would help me, but I only got garbage (-1, 0 and 1) back for the variables I would expect contain r,s,t. It worked with the element number, however.
Thanks, Markus
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Hi, what is the difference between "index" and "value"? Essentially, I am looking for the information postx provides when figuring out history points. I browsed the postx source code, located the routine and implemented it in .usr, but I expect trouble when moving to bigger domains. I also figured there must be something in nek. So I took the new hpts routine (see attached my modified .usr file), what I am getting for r,s,t for a specific point is: 1.00000000000000 -5.551115123125783E-017 -1.00000000000000 The real deal for the same point (from postx) is 4 1 1. What am I missing? Is there a way to go from the output I have to the r,s,t indices? Thanks for any help, Markus [email protected] wrote:
Markus,
those r,s,t values (not indices) are valid - do they not correspond to your expectations? Off-hand one might expect that the associated indices would be 1,NY1/2,NZ1
Paul
On Fri, 9 Jul 2010, [email protected] wrote:
Howdy Neks,
is there a function in Nek that for a given point in x,y,z would return processor number, element number and r,s,t index of the nearest GLL grid point? I was desperately trying to figure out how the interpolation routines (f.e. for history points) would help me, but I only got garbage (-1, 0 and 1) back for the variables I would expect contain r,s,t. It worked with the element number, however.
Thanks, Markus
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
C----------------------------------------------------------------------- C nek5000 user-file template C C user specified routines: C - userbc : boundary conditions C - useric : initial conditions C - uservp : variable properties C - userf : local acceleration term for fluid C - userq : local source term for scalars C - userchk: general purpose routine for checking errors etc. C C----------------------------------------------------------------------- subroutine uservp(ix,iy,iz,eg) ! set variable properties include 'SIZE' include 'TOTAL' include 'NEKUSE' integer e,f,eg c e = gllel(eg) udiff = 0.0 utrans = 0.0 return end c----------------------------------------------------------------------- subroutine userf(ix,iy,iz,eg) ! set acceleration term c c Note: this is an acceleration term, NOT a force! c Thus, ffx will subsequently be multiplied by rho(x,t). c include 'SIZE' include 'TOTAL' include 'NEKUSE' integer e,f,eg c e = gllel(eg) ffx = 0.0 ffy = 0.0 ffz = 0.0 return end c----------------------------------------------------------------------- subroutine userq(ix,iy,iz,eg) ! set source term include 'SIZE' include 'TOTAL' include 'NEKUSE' integer e,f,eg c e = gllel(eg) qvol = 0.0 source = 0.0 return end c----------------------------------------------------------------------- subroutine userbc(ix,iy,iz,iside,ieg) ! set up boundary conditions include 'SIZE' include 'TOTAL' include 'NEKUSE' c if (cbc(iside,gllel(ieg),ifield).eq.'v01') ux = 0.0 uy = 0.0 uz = 0.0 temp = 0.0 return end c----------------------------------------------------------------------- subroutine useric(ix,iy,iz,ieg) ! set up initial conditions include 'SIZE' include 'TOTAL' include 'NEKUSE' ux = 0.0 uy = 1.0 uz = 0.0 temp = 0.0 return end c----------------------------------------------------------------------- subroutine userchk() include 'SIZE' include 'TOTAL' call hptsnew() call exitt return end c----------------------------------------------------------------------- subroutine usrdat() ! This routine to modify element vertices include 'SIZE' include 'TOTAL' return end c----------------------------------------------------------------------- subroutine usrdat2() ! This routine to modify mesh coordinates include 'SIZE' include 'TOTAL' return end c----------------------------------------------------------------------- subroutine usrdat3() include 'SIZE' include 'TOTAL' return end c----------------------------------------------------------------------- subroutine hptsnew c INCLUDE 'SIZE' INCLUDE 'TOTAL' real pts(ldim,lhis) real fieldout(ldim+ldimt+1,lhis) real dist(lhis) real rst(lhis*ldim) integer rcode(lhis),elid(lhis),proc(lhis) common /hpts_r/ rst common /hpts_i/ rcode,elid,proc common /scrcg/ pm1 (lx1,ly1,lz1,lelv) ! mapped pressure common /outtmp/ wrk(lx1*ly1*lz1*lelt,ldim+ldimt+1) logical iffind integer icalld,npoints save icalld,npoints data icalld /0/ data npoints /0/ save inth_hpts nxyz = nx1*ny1*nz1 ntot = nxyz*nelt if(nid.eq.0) write(6,*) 'dump history points' if(nelgt.ne.nelgv) then if(nid.eq.0) write(6,*) & 'ABORT: hpts() no support for nelgt.ne.nelgv!' call exitt endif if(icalld.eq.0) then npoints=2 pts(1,1)=1.0 pts(2,1)=0.5 pts(3,1)=0.0 pts(1,2)=-10.0 pts(2,2)=1.0 pts(3,2)=3.0 call intpts_setup(-1.0,inth_hpts) ! use default tolerance endif if(npoints.gt.lhis) then if(nid.eq.0) write(6,*) & 'ABORT: lhis too low, increase in SIZE', npoints, lhis call exitt endif ! interpolate if(icalld.eq.0) then call findpts(inth_hpts,rcode,1, & proc,1, & elid,1, & rst,ndim, & dist,1, & pts(1,1),ndim, & pts(2,1),ndim, & pts(3,1),ndim,npoints) print *,rst(1),rst(2),rst(3) do i=1,npoints ! check return code if(rcode(i).eq.1) then if(dist(i).gt.1e-12) then write(6,'(A,4E15.7)') & ' WARNING: point on boundary or outside the mesh xy[z]d^2:' & ,(pts(k,i),k=1,ndim),dist(i) endif elseif(rcode(i).eq.2) then nfail = nfail + 1 write(6,'(A,3E15.7)') & ' WARNING: point not within mesh xy[z]: !', & (pts(k,i),k=1,ndim) endif enddo icalld = 1 endif return end c c automatically added by makenek subroutine usrsetvert(glo_num,nel,nx,ny,nz) ! to modify glo_num integer*8 glo_num(1) return end
Hi Markus, Functions in the spectral element method are expressed as: e e [1] u (r,s,t) = sum u h (r) h (s) h (t) ijk ijk i j k Where u^e_ijk are the basis coefficients for element e, with i,j,k ranging from 0 to N (1 to lx1 in the code), and the h_i's are polynomials of degree N on the interval from [-1,1]. Given r,s,t and the set of basis coefficients, one can compute u^e(r,s,t) using [1]. Note that the h_i's are Lagrangian interpolants on the Gauss-Lobatto-Legendre points, {xi_i}, i=0,...,N, such that h_i (xi_j) = delta_ij, ( = 1 for i=j, 0 otherwise ). Thus for any point (r,s,t) = (xi_i,xi_j,xi_k), we can associate basis coefficient u^e_ijk with u^e at that point. More specifically, we can associated the index set (i,j,k) with that point. However, if the continuous variables (r,s,t) do not coincide with the quadrature points then there is no associated index --- that is, the value of u^e(r,s,t) would then depend on _all_ the basis coefficients, using Eq. [1] above. The interpolation routine returns the continuous variables - as are needed for the general interpolation formula [1]. To get the indices, you could check the rst values against the GLL coordinate points that are stored in WZ in the array zgm1(), for instance. Does this help ? Paul On Fri, 9 Jul 2010, [email protected] wrote:
Hi,
what is the difference between "index" and "value"? Essentially, I am looking for the information postx provides when figuring out history points. I browsed the postx source code, located the routine and implemented it in .usr, but I expect trouble when moving to bigger domains. I also figured there must be something in nek.
So I took the new hpts routine (see attached my modified .usr file), what I am getting for r,s,t for a specific point is: 1.00000000000000 -5.551115123125783E-017 -1.00000000000000
The real deal for the same point (from postx) is 4 1 1.
What am I missing? Is there a way to go from the output I have to the r,s,t indices?
Thanks for any help, Markus
[email protected] wrote:
Markus,
those r,s,t values (not indices) are valid - do they not correspond to your expectations? Off-hand one might expect that the associated indices would be 1,NY1/2,NZ1
Paul
On Fri, 9 Jul 2010, [email protected] wrote:
Howdy Neks,
is there a function in Nek that for a given point in x,y,z would return processor number, element number and r,s,t index of the nearest GLL grid point? I was desperately trying to figure out how the interpolation routines (f.e. for history points) would help me, but I only got garbage (-1, 0 and 1) back for the variables I would expect contain r,s,t. It worked with the element number, however.
Thanks, Markus
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Hi, yes, that helped, thank you. That also explains why I saw the original point coordinates when I replaced vx, vy&vz with xm1,ym1,zm1 in the interpolation routine - the hpts function does not write out values at the nearest grid point (as the old one did), but instead interpolates all field values at the desired, arbitrary points. I think what I'll do for the sake of simplicity is to find the element number from findpts and then use the routine from postx (loops over lx1*ly1*lz1 for the specified element) to find the gll indices. These values are needed only once and thus I don't have to pay too close attention to speed. Thanks again, Markus [email protected] wrote:
Hi Markus,
Functions in the spectral element method are expressed as:
e e [1] u (r,s,t) = sum u h (r) h (s) h (t) ijk ijk i j k
Where u^e_ijk are the basis coefficients for element e, with i,j,k ranging from 0 to N (1 to lx1 in the code), and the h_i's are polynomials of degree N on the interval from [-1,1]. Given r,s,t and the set of basis coefficients, one can compute u^e(r,s,t) using [1].
Note that the h_i's are Lagrangian interpolants on the Gauss-Lobatto-Legendre points, {xi_i}, i=0,...,N, such that
h_i (xi_j) = delta_ij, ( = 1 for i=j, 0 otherwise ).
Thus for any point (r,s,t) = (xi_i,xi_j,xi_k), we can associate basis coefficient u^e_ijk with u^e at that point. More specifically, we can associated the index set (i,j,k) with that point.
However, if the continuous variables (r,s,t) do not coincide with the quadrature points then there is no associated index --- that is, the value of u^e(r,s,t) would then depend on _all_ the basis coefficients, using Eq. [1] above.
The interpolation routine returns the continuous variables - as are needed for the general interpolation formula [1].
To get the indices, you could check the rst values against the GLL coordinate points that are stored in WZ in the array zgm1(), for instance.
Does this help ?
Paul
On Fri, 9 Jul 2010, [email protected] wrote:
Hi,
what is the difference between "index" and "value"? Essentially, I am looking for the information postx provides when figuring out history points. I browsed the postx source code, located the routine and implemented it in .usr, but I expect trouble when moving to bigger domains. I also figured there must be something in nek.
So I took the new hpts routine (see attached my modified .usr file), what I am getting for r,s,t for a specific point is: 1.00000000000000 -5.551115123125783E-017 -1.00000000000000
The real deal for the same point (from postx) is 4 1 1.
What am I missing? Is there a way to go from the output I have to the r,s,t indices?
Thanks for any help, Markus
[email protected] wrote:
Markus,
those r,s,t values (not indices) are valid - do they not correspond to your expectations? Off-hand one might expect that the associated indices would be 1,NY1/2,NZ1
Paul
On Fri, 9 Jul 2010, [email protected] wrote:
Howdy Neks,
is there a function in Nek that for a given point in x,y,z would return processor number, element number and r,s,t index of the nearest GLL grid point? I was desperately trying to figure out how the interpolation routines (f.e. for history points) would help me, but I only got garbage (-1, 0 and 1) back for the variables I would expect contain r,s,t. It worked with the element number, however.
Thanks, Markus
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Hi again, Regarding finding element numbers (such as what gllel would return) for arbitrary history points, I was experimenting with the history point routine hpts a little and it seems that this function call: call findpts(inth_hpts,rcode,1, & proc,1, & elid,1, & rst,ndim, & dist,1, & pts(1,1),ndim, & pts(2,1),ndim, & pts(3,1),ndim,npoints) does not return the same element number (elid) as it can be found in the .rea file, for instance. For a handful of manually set points, I checked that the returned element number does not provide the right coordinates to enclose the desired point. Is there a way to transform the function return or am I missing something else? Thanks, Markus
Hi Markus, It's possible that it's returning the local (to the processor) element number. You would get the global (to the computational domain, defined by the .rea file) element number via: integer e,eg eg = lglel(e) Paul On Sun, 11 Jul 2010, [email protected] wrote:
Hi again,
Regarding finding element numbers (such as what gllel would return) for arbitrary history points, I was experimenting with the history point routine hpts a little and it seems that this function call:
call findpts(inth_hpts,rcode,1, & proc,1, & elid,1, & rst,ndim, & dist,1, & pts(1,1),ndim, & pts(2,1),ndim, & pts(3,1),ndim,npoints)
does not return the same element number (elid) as it can be found in the .rea file, for instance. For a handful of manually set points, I checked that the returned element number does not provide the right coordinates to enclose the desired point.
Is there a way to transform the function return or am I missing something else?
Thanks, Markus
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Hi, that was not it, I'm afraid (though I thought gllel transforms from local to global). But I also realized that the hpts routine does not deliver the correct values for me. To test this, I set up a simple 3d box with inflow, outflow and symmetry all around, Re very low, initial condition=boundary condition and 4 points. Please find all files attached, the output for the first point is OK, but then the values seem to shift positions and don't make sense. Also, when running it for 1 step or more, I get the error message: " dump history points done :: dump history points hisfile: /Users/m0s1978/nektonruns/rectest/rectest.his ERROR: .sch file already exists. ierr= 1 " And lastly, the newest repo won't compile, here is the error message: " Error: /Users/m0s1978/nek5trial/trunk/nek/drive2.f, line 1551: This name has already been assigned a data type. [NVTOT] integer*8 ntot,nvtot,nptot ---------------------^ " Thanks for any help, Markus [email protected] wrote:
Hi Markus,
It's possible that it's returning the local (to the processor) element number. You would get the global (to the computational domain, defined by the .rea file) element number via:
integer e,eg
eg = lglel(e)
Paul
On Sun, 11 Jul 2010, [email protected] wrote:
Hi again,
Regarding finding element numbers (such as what gllel would return) for arbitrary history points, I was experimenting with the history point routine hpts a little and it seems that this function call:
call findpts(inth_hpts,rcode,1, & proc,1, & elid,1, & rst,ndim, & dist,1, & pts(1,1),ndim, & pts(2,1),ndim, & pts(3,1),ndim,npoints)
does not return the same element number (elid) as it can be found in the .rea file, for instance. For a handful of manually set points, I checked that the returned element number does not provide the right coordinates to enclose the desired point.
Is there a way to transform the function return or am I missing something else?
Thanks, Markus
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Hi Markus, lglel reads "local to global element number" so input is local element number and output is global element number. I'll be seeing Stefan this week so I'm sure we can straighten out the interpolation routines. You're right about the other points - the ntotv wasn't flagged by my compiler and the .his file was taken care of by my own scripts, so I didn't catch either of them. Sorry about that. I've updated the source so that those issues should be resolved. Paul On Mon, 12 Jul 2010, [email protected] wrote:
Hi,
that was not it, I'm afraid (though I thought gllel transforms from local to global).
But I also realized that the hpts routine does not deliver the correct values for me. To test this, I set up a simple 3d box with inflow, outflow and symmetry all around, Re very low, initial condition=boundary condition and 4 points. Please find all files attached, the output for the first point is OK, but then the values seem to shift positions and don't make sense.
Also, when running it for 1 step or more, I get the error message: " dump history points done :: dump history points hisfile: /Users/m0s1978/nektonruns/rectest/rectest.his
ERROR: .sch file already exists. ierr= 1 "
And lastly, the newest repo won't compile, here is the error message: " Error: /Users/m0s1978/nek5trial/trunk/nek/drive2.f, line 1551: This name has already been assigned a data type. [NVTOT] integer*8 ntot,nvtot,nptot ---------------------^ "
Thanks for any help, Markus
[email protected] wrote:
Hi Markus,
It's possible that it's returning the local (to the processor) element number. You would get the global (to the computational domain, defined by the .rea file) element number via:
integer e,eg
eg = lglel(e)
Paul
On Sun, 11 Jul 2010, [email protected] wrote:
Hi again,
Regarding finding element numbers (such as what gllel would return) for arbitrary history points, I was experimenting with the history point routine hpts a little and it seems that this function call:
call findpts(inth_hpts,rcode,1, & proc,1, & elid,1, & rst,ndim, & dist,1, & pts(1,1),ndim, & pts(2,1),ndim, & pts(3,1),ndim,npoints)
does not return the same element number (elid) as it can be found in the .rea file, for instance. For a handful of manually set points, I checked that the returned element number does not provide the right coordinates to enclose the desired point.
Is there a way to transform the function return or am I missing something else?
Thanks, Markus
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
fixed in r548. Stefan On Jul 12, 2010, at 10:49 AM, [email protected] wrote:
Hi Markus,
lglel reads "local to global element number" so input is local element number and output is global element number.
I'll be seeing Stefan this week so I'm sure we can straighten out the interpolation routines.
You're right about the other points - the ntotv wasn't flagged by my compiler and the .his file was taken care of by my own scripts, so I didn't catch either of them. Sorry about that. I've updated the source so that those issues should be resolved.
Paul
On Mon, 12 Jul 2010, [email protected] wrote:
Hi,
that was not it, I'm afraid (though I thought gllel transforms from local to global).
But I also realized that the hpts routine does not deliver the correct values for me. To test this, I set up a simple 3d box with inflow, outflow and symmetry all around, Re very low, initial condition=boundary condition and 4 points. Please find all files attached, the output for the first point is OK, but then the values seem to shift positions and don't make sense.
Also, when running it for 1 step or more, I get the error message: " dump history points done :: dump history points hisfile: /Users/m0s1978/nektonruns/rectest/rectest.his ERROR: .sch file already exists. ierr= 1 "
And lastly, the newest repo won't compile, here is the error message: " Error: /Users/m0s1978/nek5trial/trunk/nek/drive2.f, line 1551: This name has already been assigned a data type. [NVTOT] integer*8 ntot,nvtot,nptot ---------------------^ "
Thanks for any help, Markus
[email protected] wrote:
Hi Markus, It's possible that it's returning the local (to the processor) element number. You would get the global (to the computational domain, defined by the .rea file) element number via:
integer e,eg
eg = lglel(e) Paul On Sun, 11 Jul 2010, [email protected] wrote:
Hi again, Regarding finding element numbers (such as what gllel would return) for arbitrary history points, I was experimenting with the history point routine hpts a little and it seems that this function call:
call findpts(inth_hpts,rcode,1, & proc,1, & elid,1, & rst,ndim, & dist,1, & pts(1,1),ndim, & pts(2,1),ndim, & pts(3,1),ndim,npoints) does not return the same element number (elid) as it can be found in the .rea file, for instance. For a handful of manually set points, I checked that the returned element number does not provide the right coordinates to enclose the desired point. Is there a way to transform the function return or am I missing something else? Thanks, Markus _______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Hi Paul, To expand on this topic, I am wondering what the convention being used to number the GL points for an element is? Could you explain this? Thinking of it terms of the r-s-t diagram for a nekton element, I know from previous examples that face 5 contains points 1 to lx1*ly1. Also, face 6 should contain points (lz1-1)*lx1*ly1 to lx1*ly1*lz1. What I am trying to figure out, is what points are on the remaining faces, and understanding the convention I think would be useful. Thanks for the help! - Michael ----- Original Message ----- From: [email protected] To: [email protected] Sent: Monday, July 12, 2010 5:13:47 AM GMT -06:00 US/Canada Central Subject: Re: [Nek5000-users] r,s,t index for given point fixed in r548. Stefan On Jul 12, 2010, at 10:49 AM, [email protected] wrote:
Hi Markus,
lglel reads "local to global element number" so input is local element number and output is global element number.
I'll be seeing Stefan this week so I'm sure we can straighten out the interpolation routines.
You're right about the other points - the ntotv wasn't flagged by my compiler and the .his file was taken care of by my own scripts, so I didn't catch either of them. Sorry about that. I've updated the source so that those issues should be resolved.
Paul
On Mon, 12 Jul 2010, [email protected] wrote:
Hi,
that was not it, I'm afraid (though I thought gllel transforms from local to global).
But I also realized that the hpts routine does not deliver the correct values for me. To test this, I set up a simple 3d box with inflow, outflow and symmetry all around, Re very low, initial condition=boundary condition and 4 points. Please find all files attached, the output for the first point is OK, but then the values seem to shift positions and don't make sense.
Also, when running it for 1 step or more, I get the error message: " dump history points done :: dump history points hisfile: /Users/m0s1978/nektonruns/rectest/rectest.his ERROR: .sch file already exists. ierr= 1 "
And lastly, the newest repo won't compile, here is the error message: " Error: /Users/m0s1978/nek5trial/trunk/nek/drive2.f, line 1551: This name has already been assigned a data type. [NVTOT] integer*8 ntot,nvtot,nptot ---------------------^ "
Thanks for any help, Markus
[email protected] wrote:
Hi Markus, It's possible that it's returning the local (to the processor) element number. You would get the global (to the computational domain, defined by the .rea file) element number via:
integer e,eg
eg = lglel(e) Paul On Sun, 11 Jul 2010, [email protected] wrote:
Hi again, Regarding finding element numbers (such as what gllel would return) for arbitrary history points, I was experimenting with the history point routine hpts a little and it seems that this function call:
call findpts(inth_hpts,rcode,1, & proc,1, & elid,1, & rst,ndim, & dist,1, & pts(1,1),ndim, & pts(2,1),ndim, & pts(3,1),ndim,npoints) does not return the same element number (elid) as it can be found in the .rea file, for instance. For a handful of manually set points, I checked that the returned element number does not provide the right coordinates to enclose the desired point. Is there a way to transform the function return or am I missing something else? Thanks, Markus _______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Hi Michael, Points in the volume are stored in lexicographical order: real u(lx1,ly1,lz1) implies that u(i+1,j,k) follows u(i,j,k) in memory, so the data layout (in memory) is: u(1,1,1) u(2,1,1)...u(lx1,1,1) u(1,2,1) u(2,2,1)... u(lx1,ly1,lz1) For u(i,j,k) the i index is associated with r, j with s, and k with t. There are two numberings of the faces. The one that you're most likely to interact with, and which is used for assigning bcs, unit normals, surface Jacobians, etc. is what we refer to as the "preprocessor notation". In 2D, preprocessor notation uses the counter-clockwise ordering: ^ s | Face 3 +-----------+ | | | | | | Face 4 | +-----|--> r | | | | Face 2 +-----------+ Face 1 In 3D, it's the same, save that you also have Face 5 assocated with t=-1 and Face 6 associated with t=+1. When you traverse a face, the indices describing that 2D manifold advance in positive order and in that case the associated face objects (e.g., unit normals and area) will have unit stride. There are two ways to stride across a face f=1,...,6 (and there are many routines implementing these). Here is one: integer e,f call facind (i0,i1,j0,j1,k0,k1,f) ia = 0 flux = 0 do k=k0,k1 do j=j0,j1 do i=i0,i1 ia = ia + 1 flux = flux + ( vx(i,j,k,e)*unx(ia,1,f,e) $ vy(i,j,k,e)*uny(ia,1,f,e) $ vz(i,j,k,e)*unz(ia,1,f,e) )*area(ia,1,f,e) enddo enddo enddo Here is another: call facind2 (js1,jf1,jskip1,js2,jf2,jskip2,f) ia = 0 flux = 0 do j2=js2,jf2,jskip2 do j1=js1,jf1,jskip1 ia = ia+1 flux = flux + ( vx(j1,j2,1,e)*unx(ia,1,f,e) $ + vy(j1,j2,1,e)*uny(ia,1,f,e) $ + vz(j1,j2,1,e)*unz(ia,1,f,e) )*area(ia,1,f,e) enddo enddo I like the second approach because it is faster. Note that unx,uny,unz are _face_ based structures, while vx vy vz are volumetric. I think if you study these loops for a bit you should be able to see the connections between the two data layouts. Paul On Mon, 12 Jul 2010, [email protected] wrote: Hi Paul, To expand on this topic, I am wondering what the convention being used to number the GL points for an element is? Could you explain this? Thinking of it terms of the r-s-t diagram for a nekton element, I know from previous examples that face 5 contains points 1 to lx1*ly1. Also, face 6 should contain points (lz1-1)*lx1*ly1 to lx1*ly1*lz1. What I am trying to figure out, is what points are on the remaining faces, and understanding the convention I think would be useful. Thanks for the help!
- Michael
----- Original Message ----- From: [email protected] To: [email protected] Sent: Monday, July 12, 2010 5:13:47 AM GMT -06:00 US/Canada Central Subject: Re: [Nek5000-users] r,s,t index for given point
fixed in r548. Stefan
On Jul 12, 2010, at 10:49 AM, [email protected] wrote:
Hi Markus,
lglel reads "local to global element number" so input is local element number and output is global element number.
I'll be seeing Stefan this week so I'm sure we can straighten out the interpolation routines.
You're right about the other points - the ntotv wasn't flagged by my compiler and the .his file was taken care of by my own scripts, so I didn't catch either of them. Sorry about that. I've updated the source so that those issues should be resolved.
Paul
On Mon, 12 Jul 2010, [email protected] wrote:
Hi,
that was not it, I'm afraid (though I thought gllel transforms from local to global).
But I also realized that the hpts routine does not deliver the correct values for me. To test this, I set up a simple 3d box with inflow, outflow and symmetry all around, Re very low, initial condition=boundary condition and 4 points. Please find all files attached, the output for the first point is OK, but then the values seem to shift positions and don't make sense.
Also, when running it for 1 step or more, I get the error message: " dump history points done :: dump history points hisfile: /Users/m0s1978/nektonruns/rectest/rectest.his ERROR: .sch file already exists. ierr= 1 "
And lastly, the newest repo won't compile, here is the error message: " Error: /Users/m0s1978/nek5trial/trunk/nek/drive2.f, line 1551: This name has already been assigned a data type. [NVTOT] integer*8 ntot,nvtot,nptot ---------------------^ "
Thanks for any help, Markus
[email protected] wrote:
Hi Markus, It's possible that it's returning the local (to the processor) element number. You would get the global (to the computational domain, defined by the .rea file) element number via:
integer e,eg
eg = lglel(e) Paul On Sun, 11 Jul 2010, [email protected] wrote:
Hi again, Regarding finding element numbers (such as what gllel would return) for arbitrary history points, I was experimenting with the history point routine hpts a little and it seems that this function call:
call findpts(inth_hpts,rcode,1, & proc,1, & elid,1, & rst,ndim, & dist,1, & pts(1,1),ndim, & pts(2,1),ndim, & pts(3,1),ndim,npoints) does not return the same element number (elid) as it can be found in the .rea file, for instance. For a handful of manually set points, I checked that the returned element number does not provide the right coordinates to enclose the desired point. Is there a way to transform the function return or am I missing something else? Thanks, Markus _______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Hi Paul, I am looking through the routine facind, and am wondering: Does this routine provide the loop limits for a particular face? Meaning, if I want to look at the points on face 1, the values it returns are what I loop through for r, s, and t? Also, I can't seem to find the subroutine for facind2 to have a look at it. I have found many calls for it though...where is it located? - Michael ----- Original Message ----- From: [email protected] To: [email protected] Sent: Tuesday, July 13, 2010 4:38:28 AM GMT -06:00 Guadalajara / Mexico City / Monterrey Subject: Re: [Nek5000-users] r,s,t index for given point Hi Michael, Points in the volume are stored in lexicographical order: real u(lx1,ly1,lz1) implies that u(i+1,j,k) follows u(i,j,k) in memory, so the data layout (in memory) is: u(1,1,1) u(2,1,1)...u(lx1,1,1) u(1,2,1) u(2,2,1)... u(lx1,ly1,lz1) For u(i,j,k) the i index is associated with r, j with s, and k with t. There are two numberings of the faces. The one that you're most likely to interact with, and which is used for assigning bcs, unit normals, surface Jacobians, etc. is what we refer to as the "preprocessor notation". In 2D, preprocessor notation uses the counter-clockwise ordering: ^ s | Face 3 +-----------+ | | | | | | Face 4 | +-----|--> r | | | | Face 2 +-----------+ Face 1 In 3D, it's the same, save that you also have Face 5 assocated with t=-1 and Face 6 associated with t=+1. When you traverse a face, the indices describing that 2D manifold advance in positive order and in that case the associated face objects (e.g., unit normals and area) will have unit stride. There are two ways to stride across a face f=1,...,6 (and there are many routines implementing these). Here is one: integer e,f call facind (i0,i1,j0,j1,k0,k1,f) ia = 0 flux = 0 do k=k0,k1 do j=j0,j1 do i=i0,i1 ia = ia + 1 flux = flux + ( vx(i,j,k,e)*unx(ia,1,f,e) $ vy(i,j,k,e)*uny(ia,1,f,e) $ vz(i,j,k,e)*unz(ia,1,f,e) )*area(ia,1,f,e) enddo enddo enddo Here is another: call facind2 (js1,jf1,jskip1,js2,jf2,jskip2,f) ia = 0 flux = 0 do j2=js2,jf2,jskip2 do j1=js1,jf1,jskip1 ia = ia+1 flux = flux + ( vx(j1,j2,1,e)*unx(ia,1,f,e) $ + vy(j1,j2,1,e)*uny(ia,1,f,e) $ + vz(j1,j2,1,e)*unz(ia,1,f,e) )*area(ia,1,f,e) enddo enddo I like the second approach because it is faster. Note that unx,uny,unz are _face_ based structures, while vx vy vz are volumetric. I think if you study these loops for a bit you should be able to see the connections between the two data layouts. Paul On Mon, 12 Jul 2010, [email protected] wrote: Hi Paul, To expand on this topic, I am wondering what the convention being used to number the GL points for an element is? Could you explain this? Thinking of it terms of the r-s-t diagram for a nekton element, I know from previous examples that face 5 contains points 1 to lx1*ly1. Also, face 6 should contain points (lz1-1)*lx1*ly1 to lx1*ly1*lz1. What I am trying to figure out, is what points are on the remaining faces, and understanding the convention I think would be useful. Thanks for the help! > > - Michael > > > ----- Original Message ----- > From: [email protected] > To: [email protected] > Sent: Monday, July 12, 2010 5:13:47 AM GMT -06:00 US/Canada Central > Subject: Re: [Nek5000-users] r,s,t index for given point > > fixed in r548. > Stefan > > On Jul 12, 2010, at 10:49 AM, [email protected] wrote: > >> >> Hi Markus, >> >> lglel reads "local to global element number" so input is >> local element number and output is global element number. >> >> I'll be seeing Stefan this week so I'm sure we can straighten >> out the interpolation routines. >> >> You're right about the other points - the ntotv wasn't flagged >> by my compiler and the .his file was taken care of by my own >> scripts, so I didn't catch either of them. Sorry about that. >> I've updated the source so that those issues should be >> resolved. >> >> Paul >> >> >> On Mon, 12 Jul 2010, [email protected] wrote: >> >>> Hi, >>> >>> that was not it, I'm afraid (though I thought gllel transforms from local to global). >>> >>> But I also realized that the hpts routine does not deliver t he correct values for me. To test this, I set up a simple 3d box with inflow, outflow and symmetry all around, Re very low, initial condition=boundary condition and 4 points. >>> Please find all files attached, the output for the first point is OK, but then the values seem to shift positions and don't make sense. >>> >>> Also, when running it for 1 step or more, I get the error message: >>> " >>> dump history points >>> done :: dump history points >>> hisfile: >>> /Users/m0s1978/nektonruns/rectest/rectest.his >>> ERROR: .sch file already exists. ierr= 1 >>> " >>> >>> And lastly, the newest repo won't compile, here is the error message: >>> " >>> Error: /Users/m0s1978/nek5trial/trunk/nek/drive2.f, line 1551: This name has already been assigned a data type. [NVTOT] >>> integer*8 ntot,nvtot,nptot >>> ---------------------^ >>> " >>> >>> Thanks for any help, >>> Markus >>> >>> >>> [email protected] wrote: >>>> Hi Markus, >>>> It's possible that it's returning the local (to the processor) >>>> element number. You would get the global (to the computational >>>> domain, defined by the .rea file) element number via: >>>> >>>> integer e,eg >>>> >>>> eg = lglel(e) >>>> Paul >>>> On Sun, 11 Jul 2010, [email protected] wrote: >>>>> Hi again, >>>>> Regarding finding element numbers (such as what gllel would return) for arbitrary history points, I was experimenting with the history point routine hpts a little and it seems that this function call: >>>>> >>>>> call findpts(inth_hpts,rcode,1, >>>>> & proc,1, >>>>> & elid,1, >>>>> & rst,ndim, >>>>> & dist,1, >>>>> & pts(1,1),ndim, >>>>> & pts(2,1),ndim, >>>>> & pts(3,1),ndim,npoints) >>>>> does not return the same element number (elid) as it can be found in the .rea file, for instance. >>>>> For a handful of manually set points, I checked that the returned element number does not provide the right coordinates to enclose the desired point. >>>>> Is there a way to transform the function return or am I missing something else? >>>>> Thanks, >>>>> Markus >>>>> _______________________________________________ >>>>> Nek5000-users mailing list >>>>> [email protected] >>>>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users >>>> _______________________________________________ >>>> Nek5000-users mailing list >>>> [email protected] >>>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users >>> >> _______________________________________________ >> Nek5000-users mailing list >> [email protected] >> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users > > _______________________________________________ > Nek5000-users mailing list > [email protected] > https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users > _______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Hi Michael, [NOTE: I've put the prior email in Section 14 (p.33) of www.mcs.anl.gov/~fischer/primer.pdf because it appears that my earlier email may have been garbled through deletion of all carriage returns... While there, also see Sec. 4 (p. 14) ] To locate facind2 in the source, I would use: grep -i subroutine *.f | grep -i facind2 Yes, facind gives the indices for face f. You must structure the loops exactly as below if you want to match entries on face f of volume data (in this example, vx,vy,vz) to any of the nek face-based data structures (here, unx,uny,unz,area). Paul call facind (i0,i1,j0,j1,k0,k1,f) ia = 0 flux = 0 do k=k0,k1 do j=j0,j1 do i=i0,i1 ia = ia + 1 flux = flux + ( vx(i,j,k,e)*unx(ia,1,f,e) $ vy(i,j,k,e)*uny(ia,1,f,e) $ vz(i,j,k,e)*unz(ia,1,f,e) )*area(ia,1,f,e) enddo enddo enddo On Tue, 13 Jul 2010, [email protected] wrote: Hi Paul, I am looking through the routine facind, and am wondering: Does this routine provide the loop limits for a particular face? Meaning, if I want to look at the points on face 1, the values it returns are what I loop through for r, s, and t? Also, I can't seem to find the subroutine for facind2 to have a look at it. I have found many calls for it though...where is it located? - Michael uuuu Hi Michael, Points in the volume are stored in lexicographical order: real u(lx1,ly1,lz1) implies that u(i+1,j,k) follows u(i,j,k) in memory, so the data layout (in memory) is: u(1,1,1) u(2,1,1)...u(lx1,1,1) u(1,2,1) u(2,2,1)... u(lx1,ly1,lz1) For u(i,j,k) the i index is associated with r, j with s, and k with t. There are two numberings of the faces. The one that you're most likely to interact with, and which is used for assigning bcs, unit normals, surface Jacobians, etc. is what we refer to as the "preprocessor notation". In 2D, preprocessor notation uses the counter-clockwise ordering: ^ s | Face 3 +-----------+ | | | | | | Face 4 | +-----|--> r | | | | Face 2 +-----------+ Face 1 In 3D, it's the same, save that you also have Face 5 assocated with t=-1 and Face 6 associated with t=+1. When you traverse a face, the indices describing that 2D manifold advance in positive order and in that case the associated face objects (e.g., unit normals and area) will have unit stride. There are two ways to stride across a face f=1,...,6 (and there are many routines implementing these). Here is one: integer e,f call facind (i0,i1,j0,j1,k0,k1,f) ia = 0 flux = 0 do k=k0,k1 do j=j0,j1 do i=i0,i1 ia = ia + 1 flux = flux + ( vx(i,j,k,e)*unx(ia,1,f,e) $ vy(i,j,k,e)*uny(ia,1,f,e) $ vz(i,j,k,e)*unz(ia,1,f,e) )*area(ia,1,f,e) enddo enddo enddo Here is another: call facind2 (js1,jf1,jskip1,js2,jf2,jskip2,f) ia = 0 flux = 0 do j2=js2,jf2,jskip2 do j1=js1,jf1,jskip1 ia = ia+1 flux = flux + ( vx(j1,j2,1,e)*unx(ia,1,f,e) $ + vy(j1,j2,1,e)*uny(ia,1,f,e) $ + vz(j1,j2,1,e)*unz(ia,1,f,e) )*area(ia,1,f,e) enddo enddo I like the second approach because it is faster. Note that unx,uny,unz are _face_ based structures, while vx vy vz are volumetric. I think if you study these loops for a bit you should be able to see the connections between the two data layouts. Paul On Mon, 12 Jul 2010, [email protected] wrote: Hi Paul, To expand on this topic, I am wondering what the convention being used to number the GL points for an element is? Could you explain this? Thinking of it terms of the r-s-t diagram for a nekton element, I know from previous examples that face 5 contains points 1 to lx1*ly1. Also, face 6 should contain points (lz1-1)*lx1*ly1 to lx1*ly1*lz1. What I am trying to figure out, is what points are on the remaining faces, and understanding the convention I think would be useful. Thanks for the help!
- Michael
----- Original Message ----- From: [email protected] To: [email protected] Sent: Monday, July 12, 2010 5:13:47 AM GMT -06:00 US/Canada Central Subject: Re: [Nek5000-users] r,s,t index for given point
fixed in r548. Stefan
On Jul 12, 2010, at 10:49 AM, [email protected] wrote:
Hi Markus,
lglel reads "local to global element number" so input is local element number and output is global element number.
I'll be seeing Stefan this week so I'm sure we can straighten out the interpolation routines.
You're right about the other points - the ntotv wasn't flagged by my compiler and the .his file was taken care of by my own scripts, so I didn't catch either of them. Sorry about that. I've updated the source so that those issues should be resolved.
Paul
On Mon, 12 Jul 2010, [email protected] wrote:
Hi,
that was not it, I'm afraid (though I thought gllel transforms from local to global).
But I also realized that the hpts routine does not deliver the correct values for me. To test this, I set up a simple 3d box with inflow, outflow and symmetry all around, Re very low, initial condition=boundary condition and 4 points. Please find all files attached, the output for the first point is OK, but then the values seem to shift positions and don't make sense.
Also, when running it for 1 step or more, I get the error message: " dump history points done :: dump history points hisfile: /Users/m0s1978/nektonruns/rectest/rectest.his ERROR: .sch file already exists. ierr= 1 "
And lastly, the newest repo won't compile, here is the error message: " Error: /Users/m0s1978/nek5trial/trunk/nek/drive2.f, line 1551: This name has already been assigned a data type. [NVTOT] integer*8 ntot,nvtot,nptot ---------------------^ "
Thanks for any help, Markus
[email protected] wrote:
Hi Markus, It's possible that it's returning the local (to the processor) element number. You would get the global (to the computational domain, defined by the .rea file) element number via:
integer e,eg
eg = lglel(e) Paul On Sun, 11 Jul 2010, [email protected] wrote:
Hi again, Regarding finding element numbers (such as what gllel would return) for arbitrary history points, I was experimenting with the history point routine hpts a little and it seems that this function call:
call findpts(inth_hpts,rcode,1, & proc,1, & elid,1, & rst,ndim, & dist,1, & pts(1,1),ndim, & pts(2,1),ndim, & pts(3,1),ndim,npoints) does not return the same element number (elid) as it can be found in the .rea file, for instance. For a handful of manually set points, I checked that the returned element number does not provide the right coordinates to enclose the desired point. Is there a way to transform the function return or am I missing something else? Thanks, Markus _______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Howdy, I am not sure if it is of use for anybody else, but the element number as returned from findpts (vector elid), as for example used in hpts(), can be converted to a global element number (as it would appear in the rea file) by: " lglel(elid(i)+1) " I guess the jl/jl2 routines start counting the (local) element index from 0 instead of 1 - sounds simple here, but I probably pulled out a third of my scalp hair to figure it out :-) Also, the variables returned from findpts, even though this function needs to be called on every node in parallel, are only known to node0, which makes it necessary to do something like this if one really wants to know the global number: " call bcast(proc,lhis*isize) call bcast(elid,lhis*isize) call bcast(npoints,isize) !Don't do this before jl2 functions are done !Zero out globelnum do i=1,npoints globelnum(i)=0 enddo do i=1,npoints if (nid.eq.proc(i)) then globelnum(i)=lglel(elid(i)+1) endif enddo call igop(globelnum,w1,'+ ',lhis) " Markus [email protected] wrote:
Hi Markus,
lglel reads "local to global element number" so input is local element number and output is global element number.
I'll be seeing Stefan this week so I'm sure we can straighten out the interpolation routines.
You're right about the other points - the ntotv wasn't flagged by my compiler and the .his file was taken care of by my own scripts, so I didn't catch either of them. Sorry about that. I've updated the source so that those issues should be resolved.
Paul
On Mon, 12 Jul 2010, [email protected] wrote:
Hi,
that was not it, I'm afraid (though I thought gllel transforms from local to global).
But I also realized that the hpts routine does not deliver the correct values for me. To test this, I set up a simple 3d box with inflow, outflow and symmetry all around, Re very low, initial condition=boundary condition and 4 points. Please find all files attached, the output for the first point is OK, but then the values seem to shift positions and don't make sense.
Also, when running it for 1 step or more, I get the error message: " dump history points done :: dump history points hisfile: /Users/m0s1978/nektonruns/rectest/rectest.his ERROR: .sch file already exists. ierr= 1 "
And lastly, the newest repo won't compile, here is the error message: " Error: /Users/m0s1978/nek5trial/trunk/nek/drive2.f, line 1551: This name has already been assigned a data type. [NVTOT] integer*8 ntot,nvtot,nptot ---------------------^ "
Thanks for any help, Markus
[email protected] wrote:
Hi Markus,
It's possible that it's returning the local (to the processor) element number. You would get the global (to the computational domain, defined by the .rea file) element number via:
integer e,eg
eg = lglel(e)
Paul
On Sun, 11 Jul 2010, [email protected] wrote:
Hi again,
Regarding finding element numbers (such as what gllel would return) for arbitrary history points, I was experimenting with the history point routine hpts a little and it seems that this function call:
call findpts(inth_hpts,rcode,1, & proc,1, & elid,1, & rst,ndim, & dist,1, & pts(1,1),ndim, & pts(2,1),ndim, & pts(3,1),ndim,npoints)
does not return the same element number (elid) as it can be found in the .rea file, for instance. For a handful of manually set points, I checked that the returned element number does not provide the right coordinates to enclose the desired point.
Is there a way to transform the function return or am I missing something else?
Thanks, Markus
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Hi Markus, sorry that you had to struggle with all this. You're right the local element number is a 0-based index on the C-side. I think I will change the Fortran C-wrapper back to an index-1 to be compatible with Nek. Some time ago we changed lglel(iel,nid) to lglel(iel) to reduce the memory footprint. Moreover, we could not think of an application where we really need lglel(iel,nid). What are you trying to do? Stefan On Jul 14, 2010, at 9:55 PM, [email protected] wrote:
Howdy,
I am not sure if it is of use for anybody else, but the element number as returned from findpts (vector elid), as for example used in hpts(), can be converted to a global element number (as it would appear in the rea file) by: " lglel(elid(i)+1) " I guess the jl/jl2 routines start counting the (local) element index from 0 instead of 1 - sounds simple here, but I probably pulled out a third of my scalp hair to figure it out :-)
Also, the variables returned from findpts, even though this function needs to be called on every node in parallel, are only known to node0, which makes it necessary to do something like this if one really wants to know the global number: " call bcast(proc,lhis*isize) call bcast(elid,lhis*isize) call bcast(npoints,isize) !Don't do this before jl2 functions are done !Zero out globelnum do i=1,npoints globelnum(i)=0 enddo do i=1,npoints if (nid.eq.proc(i)) then globelnum(i)=lglel(elid(i)+1) endif enddo call igop(globelnum,w1,'+ ',lhis) "
Markus
[email protected] wrote:
Hi Markus, lglel reads "local to global element number" so input is local element number and output is global element number. I'll be seeing Stefan this week so I'm sure we can straighten out the interpolation routines. You're right about the other points - the ntotv wasn't flagged by my compiler and the .his file was taken care of by my own scripts, so I didn't catch either of them. Sorry about that. I've updated the source so that those issues should be resolved. Paul On Mon, 12 Jul 2010, [email protected] wrote:
Hi,
that was not it, I'm afraid (though I thought gllel transforms from local to global).
But I also realized that the hpts routine does not deliver the correct values for me. To test this, I set up a simple 3d box with inflow, outflow and symmetry all around, Re very low, initial condition=boundary condition and 4 points. Please find all files attached, the output for the first point is OK, but then the values seem to shift positions and don't make sense.
Also, when running it for 1 step or more, I get the error message: " dump history points done :: dump history points hisfile: /Users/m0s1978/nektonruns/rectest/rectest.his ERROR: .sch file already exists. ierr= 1 "
And lastly, the newest repo won't compile, here is the error message: " Error: /Users/m0s1978/nek5trial/trunk/nek/drive2.f, line 1551: This name has already been assigned a data type. [NVTOT] integer*8 ntot,nvtot,nptot ---------------------^ "
Thanks for any help, Markus
[email protected] wrote:
Hi Markus,
It's possible that it's returning the local (to the processor) element number. You would get the global (to the computational domain, defined by the .rea file) element number via:
integer e,eg
eg = lglel(e)
Paul
On Sun, 11 Jul 2010, [email protected] wrote:
Hi again,
Regarding finding element numbers (such as what gllel would return) for arbitrary history points, I was experimenting with the history point routine hpts a little and it seems that this function call:
call findpts(inth_hpts,rcode,1, & proc,1, & elid,1, & rst,ndim, & dist,1, & pts(1,1),ndim, & pts(2,1),ndim, & pts(3,1),ndim,npoints)
does not return the same element number (elid) as it can be found in the .rea file, for instance. For a handful of manually set points, I checked that the returned element number does not provide the right coordinates to enclose the desired point.
Is there a way to transform the function return or am I missing something else?
Thanks, Markus
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Hi Stefan, please let me know once you do change the wrapper. I think the way I am sailing around not having lglel(iel,nid) is fine since it is only done once, at the beginning of the simulation. What I was working on is a new way to do the recycling boundary condition based on coordinates instead of element counts. This is necessary for us due to the usage of Gambit for our more complex geometries. We are still in the testing phase, but I think I got it to work. What the functions I wrote do is to search for the correct recycling point for every grid point on a "v" face (for which I need findpts), then to fill the ptr/glo_num variable accordingly (thanks to Paul in the right "memory" order - that would have been another third of my hair because intuitively, I would have thought its u(1,1,2,1) after u(1,1,1,1), I guess - so thanks). Markus [email protected] wrote:
Hi Markus,
sorry that you had to struggle with all this. You're right the local element number is a 0-based index on the C-side. I think I will change the Fortran C-wrapper back to an index-1 to be compatible with Nek.
Some time ago we changed lglel(iel,nid) to lglel(iel) to reduce the memory footprint. Moreover, we could not think of an application where we really need lglel(iel,nid).
What are you trying to do?
Stefan
On Jul 14, 2010, at 9:55 PM, [email protected] wrote:
Howdy,
I am not sure if it is of use for anybody else, but the element number as returned from findpts (vector elid), as for example used in hpts(), can be converted to a global element number (as it would appear in the rea file) by: " lglel(elid(i)+1) " I guess the jl/jl2 routines start counting the (local) element index from 0 instead of 1 - sounds simple here, but I probably pulled out a third of my scalp hair to figure it out :-)
Also, the variables returned from findpts, even though this function needs to be called on every node in parallel, are only known to node0, which makes it necessary to do something like this if one really wants to know the global number: " call bcast(proc,lhis*isize) call bcast(elid,lhis*isize) call bcast(npoints,isize) !Don't do this before jl2 functions are done !Zero out globelnum do i=1,npoints globelnum(i)=0 enddo do i=1,npoints if (nid.eq.proc(i)) then globelnum(i)=lglel(elid(i)+1) endif enddo call igop(globelnum,w1,'+ ',lhis) "
Markus
[email protected] wrote:
Hi Markus, lglel reads "local to global element number" so input is local element number and output is global element number. I'll be seeing Stefan this week so I'm sure we can straighten out the interpolation routines. You're right about the other points - the ntotv wasn't flagged by my compiler and the .his file was taken care of by my own scripts, so I didn't catch either of them. Sorry about that. I've updated the source so that those issues should be resolved. Paul On Mon, 12 Jul 2010, [email protected] wrote:
Hi,
that was not it, I'm afraid (though I thought gllel transforms from local to global).
But I also realized that the hpts routine does not deliver the correct values for me. To test this, I set up a simple 3d box with inflow, outflow and symmetry all around, Re very low, initial condition=boundary condition and 4 points. Please find all files attached, the output for the first point is OK, but then the values seem to shift positions and don't make sense.
Also, when running it for 1 step or more, I get the error message: " dump history points done :: dump history points hisfile: /Users/m0s1978/nektonruns/rectest/rectest.his ERROR: .sch file already exists. ierr= 1 "
And lastly, the newest repo won't compile, here is the error message: " Error: /Users/m0s1978/nek5trial/trunk/nek/drive2.f, line 1551: This name has already been assigned a data type. [NVTOT] integer*8 ntot,nvtot,nptot ---------------------^ "
Thanks for any help, Markus
[email protected] wrote:
Hi Markus,
It's possible that it's returning the local (to the processor) element number. You would get the global (to the computational domain, defined by the .rea file) element number via:
integer e,eg
eg = lglel(e)
Paul
On Sun, 11 Jul 2010, [email protected] wrote:
Hi again,
Regarding finding element numbers (such as what gllel would return) for arbitrary history points, I was experimenting with the history point routine hpts a little and it seems that this function call:
call findpts(inth_hpts,rcode,1, & proc,1, & elid,1, & rst,ndim, & dist,1, & pts(1,1),ndim, & pts(2,1),ndim, & pts(3,1),ndim,npoints)
does not return the same element number (elid) as it can be found in the .rea file, for instance. For a handful of manually set points, I checked that the returned element number does not provide the right coordinates to enclose the desired point.
Is there a way to transform the function return or am I missing something else?
Thanks, Markus
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
HI Markus -- give me a few minutes and I will send out an email with a test case that can help you with RBC Best, Aleks On Fri, 16 Jul 2010, [email protected] wrote:
Hi Stefan,
please let me know once you do change the wrapper.
I think the way I am sailing around not having lglel(iel,nid) is fine since it is only done once, at the beginning of the simulation.
What I was working on is a new way to do the recycling boundary condition based on coordinates instead of element counts. This is necessary for us due to the usage of Gambit for our more complex geometries. We are still in the testing phase, but I think I got it to work.
What the functions I wrote do is to search for the correct recycling point for every grid point on a "v" face (for which I need findpts), then to fill the ptr/glo_num variable accordingly (thanks to Paul in the right "memory" order - that would have been another third of my hair because intuitively, I would have thought its u(1,1,2,1) after u(1,1,1,1), I guess - so thanks).
Markus
[email protected] wrote:
Hi Markus,
sorry that you had to struggle with all this. You're right the local element number is a 0-based index on the C-side. I think I will change the Fortran C-wrapper back to an index-1 to be compatible with Nek. Some time ago we changed lglel(iel,nid) to lglel(iel) to reduce the memory footprint. Moreover, we could not think of an application where we really need lglel(iel,nid).
What are you trying to do?
Stefan
On Jul 14, 2010, at 9:55 PM, [email protected] wrote:
Howdy,
I am not sure if it is of use for anybody else, but the element number as returned from findpts (vector elid), as for example used in hpts(), can be converted to a global element number (as it would appear in the rea file) by: " lglel(elid(i)+1) " I guess the jl/jl2 routines start counting the (local) element index from 0 instead of 1 - sounds simple here, but I probably pulled out a third of my scalp hair to figure it out :-)
Also, the variables returned from findpts, even though this function needs to be called on every node in parallel, are only known to node0, which makes it necessary to do something like this if one really wants to know the global number: " call bcast(proc,lhis*isize) call bcast(elid,lhis*isize) call bcast(npoints,isize) !Don't do this before jl2 functions are done !Zero out globelnum do i=1,npoints globelnum(i)=0 enddo do i=1,npoints if (nid.eq.proc(i)) then globelnum(i)=lglel(elid(i)+1) endif enddo call igop(globelnum,w1,'+ ',lhis) "
Markus
[email protected] wrote:
Hi Markus, lglel reads "local to global element number" so input is local element number and output is global element number. I'll be seeing Stefan this week so I'm sure we can straighten out the interpolation routines. You're right about the other points - the ntotv wasn't flagged by my compiler and the .his file was taken care of by my own scripts, so I didn't catch either of them. Sorry about that. I've updated the source so that those issues should be resolved. Paul On Mon, 12 Jul 2010, [email protected] wrote:
Hi,
that was not it, I'm afraid (though I thought gllel transforms from local to global).
But I also realized that the hpts routine does not deliver the correct values for me. To test this, I set up a simple 3d box with inflow, outflow and symmetry all around, Re very low, initial condition=boundary condition and 4 points. Please find all files attached, the output for the first point is OK, but then the values seem to shift positions and don't make sense.
Also, when running it for 1 step or more, I get the error message: " dump history points done :: dump history points hisfile: /Users/m0s1978/nektonruns/rectest/rectest.his ERROR: .sch file already exists. ierr= 1 "
And lastly, the newest repo won't compile, here is the error message: " Error: /Users/m0s1978/nek5trial/trunk/nek/drive2.f, line 1551: This name has already been assigned a data type. [NVTOT] integer*8 ntot,nvtot,nptot ---------------------^ "
Thanks for any help, Markus
[email protected] wrote:
Hi Markus,
It's possible that it's returning the local (to the processor) element number. You would get the global (to the computational domain, defined by the .rea file) element number via:
integer e,eg
eg = lglel(e)
Paul
On Sun, 11 Jul 2010, [email protected] wrote:
Hi again,
Regarding finding element numbers (such as what gllel would return) for arbitrary history points, I was experimenting with the history point routine hpts a little and it seems that this function call:
call findpts(inth_hpts,rcode,1, & proc,1, & elid,1, & rst,ndim, & dist,1, & pts(1,1),ndim, & pts(2,1),ndim, & pts(3,1),ndim,npoints)
does not return the same element number (elid) as it can be found in the .rea file, for instance. For a handful of manually set points, I checked that the returned element number does not provide the right coordinates to enclose the desired point.
Is there a way to transform the function return or am I missing something else?
Thanks, Markus
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
_______________________________________________ Nek5000-users mailing list [email protected] https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
participants (1)
-
nek5000-users＠lists.mcs.anl.gov