[Nek5000-users] opdiv and multd

[nek5000-users at lists.mcs.anl.gov]

Hi Francesco,

Several of the routines internal to nek are the so-called weak
derivatives, implying that they somehow involve the mass matrix
or the transpose of the operator, so it takes a bit of investigation
to understand precisely what each operator is doing.  (Note that the
functionality also differs slightly between PN-PN and PN-PN-2 
formulations.)

As an example, axhelm does not multiply by the (negative) Laplacian.
Instead, it creates (for h1==1, h2==0):

    w = A * u

      = grad^T B grad u

where B is the mass matrix.

Similarly, opdiv is most likely returning   d = B*grad u (or perhaps
the negative of this).     Moreover, d is returned on the pressure
mesh (mesh 2), which is the same as the velocity mesh (mesh 1) for 
the PN-PN formulation that you use.

In Nek, the mass matrix is diagonal and the mesh 1 inverse is stored
in binvm1.  It seems probable that you'll get the correct result
after

       call col2(d,binvm1,n)

where d holds the output of opdiv and n is the number of points on
mesh 1, assuming you're using PN-Pn.

Paul



On Wed, 8 Dec 2010, nek5000-users at lists.mcs.anl.gov wrote:

> Hi all,
>
> I am having some problems using the subroutines opdiv and multd.
> For some reasons they do not give me correct results.
> Do they need any special treatment or precautions?
>
> I was trying to compute the divergence of a 3D field, and i realized that in 
> my case the multd subroutine used
> by opdiv does not gives good results, while the same derivative computed from 
> gram1 gives no problems.
>
> For example, I overwrite an analytical field :
>
>         ...
>         parameter(lt=lx1*ly1*lz1*lelt)
>         real       Res(lt)
>         ...
>                    vx(i,j,k,e)=   1.0*sin(x)*cos(y)*cos(z)
>         ...
>             call print_line(vx)
>         ...
>         iflg = 1
>         call multd (RES(1),vx,rym2,sym2,tym2,2,iflg)
>            call print_line(Res(1) )
>
>        call gradm1(work1,Res(1),work2,vx)
>            call print_line(Res(1) )
>         ...
>         if (lx2.ne.lx1) then
>         ...
>         STOP
>         endif
>
>
> and  print output a (xr,:,zr) the result  (performed by print_line) .
> In vxi.eps i compare the result (lines) with the analytical solution( 
> symbols). The result of multd is not correct while all
> the others match perfectly the analytical solution.
> In vxy.eps i plot only the result of multd. Note that on y i have 12 elements 
> the same number of peaks of the solution.
>
> Does anybody have an idea on how to fix this?
>
> Thank you very much for any help or suggestion
>
> francesco
>
>
>
>