I see a problem with the use of diffover: at the boundary it takes no
difference at all. (I also see a problem with the name, as it is not a
reduction as the other 'over' methods). Thus, I would use
diff2. Furthermore, both functions use a forward difference. One could
conceive a backward difference also, or their average, or, for data on
a uniform grid, the centered difference. Assuming periodic boundary
conditions, the centered difference could be obtained by rotating. If
the boundary conditions are not periodic then the boundary values
could be thrown away. My answer would then be:
use v5.36;
use PDL;
use PDL::NiceSlice;
sub partial($var, $f){
my $g=$f->mv($var,0);
return (($g->rotate(-1)-$g->rotate(1))/2)->mv(0, $var);
}
sub curl($f){
return
pdl(
partial(1,$f((2)))-partial(2,$f((1))),
partial(2,$f((0)))-partial(0,$f((2))),
partial(0,$f((1)))-partial(1,$f((0))),
)->mv(-1,0);
}
sub div($f){
return partial(0,$f((0)))+partial(1,$f((1)))+partial(2,$f((2)));
}
sub trim3d($f){
return $f->(:,1:-2,1:-2,1:-2);
}
For example:
my $z=zeroes(5,5,5);
my $v=pdl(-$z->yvals, $z->xvals, $z->zvals)->mv(-1,0);
say trim3d(curl($v));
yields
[
[
[
[0 0 2]
[0 0 2]
[0 0 2]
]
[
[0 0 2]
[0 0 2]
[0 0 2]
]
[
[0 0 2]
[0 0 2]
[0 0 2]
]
]
[
[
[0 0 2]
[0 0 2]
[0 0 2]
]
[
[0 0 2]
[0 0 2]
[0 0 2]
]
[
[0 0 2]
[0 0 2]
[0 0 2]
]
]
[
[
[0 0 2]
[0 0 2]
[0 0 2]
]
[
[0 0 2]
[0 0 2]
[0 0 2]
]
[
[0 0 2]
[0 0 2]
[0 0 2]
]
]
]
as curl(-y,x,z)=(0,0,2).
Regards,
Luis
On Wed, Nov 20, 2024 at 10:17:36AM +0000, Ed . wrote:
> This (untested) should work, as it is a fairly direct translation of the
> formula:
>
> $pP = $vec->slice('(0)')->diffover; # no mv as x dim already bottom
> $px = $coords->slice('(0)')->diffover;
> $pQ = $vec->slice('(1)')->mv(1,0)->diffover->mv(0,1);
> $py = $coords->slice('(1)')->mv(1,0)->diffover->mv(0,1);
> $pR = $vec->slice('(2)')->mv(2,0)->diffover->mv(0,2);
> $pz = $coords->slice('(2)')->mv(2,0)->diffover->mv(0,2);
>
> $curl = $vec->zeroes;
> $curl->slice('(0)') .= $pR/$py - $pQ/$pz;
> $curl->slice('(1)') .= $pP/$pz - $pR/$px;
> $curl->slice('(2)') .= $pQ/$px - $pP/$py;
>
> The $curl could be a cat of those 3 expressions, with ->mv(-1,0); there's
> probably a clever way to make one copy of each ndarray and then do inplace
> operations with less mving, and the whole thing could become a PP operation,
> but let's see if this is conceptually correct first!
>
> Best regards,
> Ed
>
> ________________________________
> From: Ed . <ej...@hotmail.com>
> Sent: 20 November 2024 9:54 AM
> To: perldl <pdl-gene...@lists.sourceforge.net>; pdl-devel
> <pdl-devel@lists.sourceforge.net>; Ingo Schmid <ingo...@gmx.at>
> Subject: Re: [Pdl-devel] curl of vector
>
> Hi Ingo,
>
> I'm not aware of any, but I had a quick google to find the formula (I know
> vaguely what divergence and curl are having watched a 3blue1brown video about
> it some time ago). I couldn't find any implementations in Python or Fortran.
>
> This
> (https://openstax.org/books/calculus-volume-3/pages/6-5-divergence-and-curl
> formula 6.17) indicates the formula for curl is:
>
> If F=⟨P,Q,R⟩ is a vector field in R3, and Px, Py, Pz, Qy, Qx, Qz, Rz, Rx,
> and Ry all exist, then the curl of F is defined by
>
> curl F = (Ry−Qz)i + (Pz−Rx)j + (Qx−Py)k
> = (∂R/∂y − ∂Q/∂z)i + (∂P/∂z − ∂R/∂x)j + (∂Q/∂x − ∂P/∂y)k
>
> This suggests to me that for a 3D problem, you need a coordinates ndarray dim
> (3,x,y,z) and a vector field ndarray with the same dims. You can do the
> partial differentiation numerically by some mv-ing and then diffover
> (https://metacpan.org/pod/PDL::Ufunc#diffover), then the notional i,j,k above
> obviously indicate the final components of curl vector at each point. More to
> follow after I've figured out how to do the partial stuff!
>
> Best regards,
> Ed
>
> ________________________________
> From: Ingo Schmid via pdl-devel <pdl-devel@lists.sourceforge.net>
> Sent: 19 November 2024 6:06 PM
> To: perldl <pdl-gene...@lists.sourceforge.net>; pdl-devel
> <pdl-devel@lists.sourceforge.net>
> Subject: [Pdl-devel] curl of vector
>
>
> Hi,
>
> is there any implementations of calculating the curl of a vector field around?
>
> Best wishes
>
> Ingo
> _______________________________________________
> pdl-devel mailing list
> pdl-devel@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/pdl-devel
--
o
W. Luis Mochán, | tel:(52)(777)329-1734 /<(*)
Instituto de Ciencias Físicas, UNAM | fax:(52)(777)317-5388 `>/ /\
Av. Universidad s/n CP 62210 | (*)/\/ \
Cuernavaca, Morelos, México | moc...@fis.unam.mx /\_/\__/
GPG: 791EB9EB, C949 3F81 6D9B 1191 9A16 C2DF 5F0A C52B 791E B9EB
_______________________________________________
pdl-devel mailing list
pdl-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/pdl-devel
