Howdy Nek's,<br><br>I have a question regarding axhelm (and other routines likes ophinv).<br><br>Paul wrote:<br><br><div style="margin-left:40px"><i>"As an example, axhelm does not multiply by the (negative) <span class="il">Laplacian</span>.<br>
Instead, it creates (for h1==1, h2==0):<br></i>
<i><br>
w = A * u<br></i>
<i><br>
= grad^T B grad u<br></i>
<i><br>
where B is the mass matrix."<br></i></div><br>So basically if (h1==1,h2==0), you only need to multiply by the inverse of the mass matrix to get laplacian(u) (+ some direct stiffness sumation). Am I right?<br>What about if h1==0 and h2==1? What does axhelm return in this case?<br>
<br>Same question about ophinv:<br><br><div style="margin-left:40px">subroutine ophinv (out1,out2,out3,inp1,inp2,inp3,h1,h2,tolh,nmxi)<br>C----------------------------------------------------------------------<br>C<br>C OUT = (H1*A+H2*B)-1 * INP (implicit)<br>
</div><br>I presume that if (h1==1,h2==0), it solves: laplacian(out) = input. Is that right? What about (h1==0,h2==1) again?<br><br>Regards and happy chinesse new year!<br><br><div class="gmail_quote">2010/12/14 <span dir="ltr"><<a href="mailto:nek5000-users@lists.mcs.anl.gov" target="_blank">nek5000-users@lists.mcs.anl.gov</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">dssum() is just a wrapper for call gs_op() which takes only the input field from the dssum call. But I agree in general it's a good idea to call all subroutines with the correct arguments.<br>
<br>
-Stefan<br>
<div class="HOEnZb"><div class="h5"><br>
<br>
<br>
<br>
<br>
----- Original Message -----<br>
From: <a href="mailto:nek5000-users@lists.mcs.anl.gov">nek5000-users@lists.mcs.anl.gov</a><br>
To: <a href="mailto:nek5000-users@lists.mcs.anl.gov">nek5000-users@lists.mcs.anl.gov</a><br>
Sent: Tuesday, December 14, 2010 6:50:15 PM<br>
Subject: Re: [Nek5000-users] opdiv and multd<br>
<br>
Sorry i have to do a correction.<br>
<br>
A user point it out that the subroutine is defined as dssum(u,nx,ny,nz)<br>
in dssum.f .<br>
Thus it should be called with 4 arguments :<br>
<br>
call dssum(RES(1),nx1,ny1,nz1)<br>
<br>
This did not give me problems and it seems to me that nx,ny,nz are not<br>
used.<br>
But it is better to be aware of it since i do not know how the compiler<br>
will read the wrong call.<br>
<br>
Sorry about that.<br>
francesco<br>
<br>
<br>
<br>
<br>
On 12/08/2010 04:02 PM, <a href="mailto:nek5000-users@lists.mcs.anl.gov">nek5000-users@lists.mcs.anl.gov</a> wrote:<br>
><br>
> Thanks Francesco!<br>
><br>
> Paul<br>
><br>
><br>
> On Wed, 8 Dec 2010, <a href="mailto:nek5000-users@lists.mcs.anl.gov">nek5000-users@lists.mcs.anl.gov</a> wrote:<br>
><br>
>> Thank you Paul,<br>
>> that solved the problem.<br>
>><br>
>> The correct result is obtained by :<br>
>> call multd (RES(1),vx,rym2,sym2,tym2,2,iflg)<br>
>> call dssum(RES(1))<br>
>> call col2(RES(1),binvm1,lt)<br>
>><br>
>> francesco<br>
>><br>
>><br>
>><br>
>> On 12/08/2010 02:17 PM, <a href="mailto:nek5000-users@lists.mcs.anl.gov">nek5000-users@lists.mcs.anl.gov</a> wrote:<br>
>>><br>
>>> Hi Francesco,<br>
>>><br>
>>> Several of the routines internal to nek are the so-called weak<br>
>>> derivatives, implying that they somehow involve the mass matrix<br>
>>> or the transpose of the operator, so it takes a bit of investigation<br>
>>> to understand precisely what each operator is doing. (Note that the<br>
>>> functionality also differs slightly between PN-PN and PN-PN-2<br>
>>> formulations.)<br>
>>><br>
>>> As an example, axhelm does not multiply by the (negative) Laplacian.<br>
>>> Instead, it creates (for h1==1, h2==0):<br>
>>><br>
>>> w = A * u<br>
>>><br>
>>> = grad^T B grad u<br>
>>><br>
>>> where B is the mass matrix.<br>
>>><br>
>>> Similarly, opdiv is most likely returning d = B*grad u (or perhaps<br>
>>> the negative of this). Moreover, d is returned on the pressure<br>
>>> mesh (mesh 2), which is the same as the velocity mesh (mesh 1) for<br>
>>> the PN-PN formulation that you use.<br>
>>><br>
>>> In Nek, the mass matrix is diagonal and the mesh 1 inverse is stored<br>
>>> in binvm1. It seems probable that you'll get the correct result<br>
>>> after<br>
>>><br>
>>> call col2(d,binvm1,n)<br>
>>><br>
>>> where d holds the output of opdiv and n is the number of points on<br>
>>> mesh 1, assuming you're using PN-Pn.<br>
>>><br>
>>> Paul<br>
>>><br>
>>><br>
>>><br>
>>> On Wed, 8 Dec 2010, <a href="mailto:nek5000-users@lists.mcs.anl.gov">nek5000-users@lists.mcs.anl.gov</a> wrote:<br>
>>><br>
>>>> Hi all,<br>
>>>><br>
>>>> I am having some problems using the subroutines opdiv and multd.<br>
>>>> For some reasons they do not give me correct results.<br>
>>>> Do they need any special treatment or precautions?<br>
>>>><br>
>>>> I was trying to compute the divergence of a 3D field, and i<br>
>>>> realized that in my case the multd subroutine used<br>
>>>> by opdiv does not gives good results, while the same derivative<br>
>>>> computed from gram1 gives no problems.<br>
>>>><br>
>>>> For example, I overwrite an analytical field :<br>
>>>><br>
>>>> ...<br>
>>>> parameter(lt=lx1*ly1*lz1*lelt)<br>
>>>> real Res(lt)<br>
>>>> ...<br>
>>>> vx(i,j,k,e)= 1.0*sin(x)*cos(y)*cos(z)<br>
>>>> ...<br>
>>>> call print_line(vx)<br>
>>>> ...<br>
>>>> iflg = 1<br>
>>>> call multd (RES(1),vx,rym2,sym2,tym2,2,iflg)<br>
>>>> call print_line(Res(1) )<br>
>>>><br>
>>>> call gradm1(work1,Res(1),work2,vx)<br>
>>>> call print_line(Res(1) )<br>
>>>> ...<br>
>>>> if (lx2.ne.lx1) then<br>
>>>> ...<br>
>>>> STOP<br>
>>>> endif<br>
>>>><br>
>>>><br>
>>>> and print output a (xr,:,zr) the result (performed by print_line) .<br>
>>>> In vxi.eps i compare the result (lines) with the analytical<br>
>>>> solution( symbols). The result of multd is not correct while all<br>
>>>> the others match perfectly the analytical solution.<br>
>>>> In vxy.eps i plot only the result of multd. Note that on y i have<br>
>>>> 12 elements the same number of peaks of the solution.<br>
>>>><br>
>>>> Does anybody have an idea on how to fix this?<br>
>>>><br>
>>>> Thank you very much for any help or suggestion<br>
>>>><br>
>>>> francesco<br>
>>>><br>
>>>><br>
>>>><br>
>>>><br>
>>> _______________________________________________<br>
>>> Nek5000-users mailing list<br>
>>> <a href="mailto:Nek5000-users@lists.mcs.anl.gov">Nek5000-users@lists.mcs.anl.gov</a><br>
>>> <a href="https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users" target="_blank">https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users</a><br>
>><br>
>> _______________________________________________<br>
>> Nek5000-users mailing list<br>
>> <a href="mailto:Nek5000-users@lists.mcs.anl.gov">Nek5000-users@lists.mcs.anl.gov</a><br>
>> <a href="https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users" target="_blank">https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users</a><br>
>><br>
> _______________________________________________<br>
> Nek5000-users mailing list<br>
> <a href="mailto:Nek5000-users@lists.mcs.anl.gov">Nek5000-users@lists.mcs.anl.gov</a><br>
> <a href="https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users" target="_blank">https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users</a><br>
><br>
<br>
_______________________________________________<br>
Nek5000-users mailing list<br>
<a href="mailto:Nek5000-users@lists.mcs.anl.gov">Nek5000-users@lists.mcs.anl.gov</a><br>
<a href="https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users" target="_blank">https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users</a><br>
_______________________________________________<br>
Nek5000-users mailing list<br>
<a href="mailto:Nek5000-users@lists.mcs.anl.gov">Nek5000-users@lists.mcs.anl.gov</a><br>
<a href="https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users" target="_blank">https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users</a><br>
</div></div></blockquote></div><br><br clear="all"><br>-- <br>Jean-Christophe