Implemented ReflectiveShell as optional boundary #512
Conversation
There was a problem hiding this comment.
Hey @ehlertdo,
Thanks for sharing your implementation. In general, the tool looks good for me.
I have tested it and found some deviation when using large steps. I guess this is caused by the position where you calculate the normal. This is given at the current position which is not necessarily a position on the sphere. Would you have an idea how to fix this, or can this be neglected in regular cases?
Beside this I would ask you to add your changes to the Changes.md and maybe you can provide a simple test case as a notebook and as unit test.
|
Hi Julien, |
include/crpropa/module/Boundary.h
Outdated
| @param cen vector corresponding to the center of the sphere | ||
| @param r value corresponding to the radius of the shell | ||
| */ | ||
| ReflectiveShell(Vector3d cen, double r); |
There was a problem hiding this comment.
Code Style: use full words for variable names: center instead of cen
src/module/Boundary.cpp
Outdated
| @@ -40,6 +40,58 @@ std::string PeriodicBox::getDescription() const { | |||
| return s.str(); | |||
| } | |||
|
|
|||
|
|
|||
| ReflectiveShell::ReflectiveShell(Vector3d cen, double r) : | |||
src/module/Boundary.cpp
Outdated
| Vector3d previousPosition = c->previous.getPosition(); | ||
| // get point where trajectory intersects the shell boundary | ||
| double p_half = previousPosition.dot(currentDirection) / currentDirection.getR2(); | ||
| double q = (previousPosition.getR2() - radius*radius) / currentDirection.getR2(); |
There was a problem hiding this comment.
Code Style: spaces around operators. Use a * b instead of a*b.
This also happens below. I mark it only once here.
src/module/Boundary.cpp
Outdated
| std::string ReflectiveShell::getDescription() const { | ||
| std::stringstream ss; | ||
| ss << "Reflective sphere: " << std::endl | ||
| << " Center: " << center << std::endl |
There was a problem hiding this comment.
Is it useful to define some default unit for the printing like kpc or Mpc?
test/testBreakCondition.cpp
Outdated
|
|
||
| shell.process(&c); | ||
|
|
||
| EXPECT_NEAR(89.9965, c.current.getPosition().x, 1e-4); |
There was a problem hiding this comment.
Please leave a comment in the code on how you get to those expected values (or even calculate them with the trigonometric functions here).
ehlertdo
force-pushed
the
ReflectiveShell
branch
from
50cc7ac to
a9ab582
Compare
Hi,
I needed a reflective boundary in my 3D simulations but in the form of a sphere instead of the currently available ReflectiveBox. Maybe this feature could be of more general interest since astrophysical environments are often more spherical than they are boxy.
Note that the shell reflects in both directions, i.e. outside particles are kept out and inside particles are kept in.
Best,
Domenik
PS: I copied the distance() and normal() functions from the Sphere class in Geometry.cpp but probably there's a more elegant way to do this with inheritance.