Severin
A general comment: Your variable pio has only 7 digits of accuracy. You typed more digits, but since this is a single precision constant, it is rounded to 7 digits. I recommend using "acos(-1.0d0)" instead.
Your loops look correct. If you want to test your code, then you could initialize a grid function to zero before your loop, and set it to one or two inside the loop. This will test whether your loop executes for all points, and which if branch is taken.
There might be another function called after your loop runs. For example, interpolation from a neighbouring patch might overwrite some data. You could output a few grid point values after your loop and compare with what you see in your images to check this.
-erik
On Mon, Dec 2, 2019 at 10:43 AM Severin Frank severin.frank@uni-tuebingen.de wrote:
Dear all,
I'm calculating radial oscillations on a Neutron Star within a multipatch (Llama) environment. The coordinate system is Thornburg04. For my time evolution of the fluid element displacement equations, I'm scheduling routines for my boundary conditions in MoL_PostStep. The conditions are with respect to the (changing) surface of the star TOV_surface. I want to loop over all grid points to apply the conditions on certain grid points (close to the surface) in Fortran.
!$OMP PARALLEL DO private(i,j,k) do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) if (grid_r(i,j,k) >= (TOV_surface-rprec) .AND. grid_r(i,j,k) <= (TOV_surface+ rprec)) then drXi(i,j,k) = 0 Pidot(i,j,k) = Wrinv(i,j,k) * P(i,j,k) * drrXi(i,j,k) + Qr(i,j,k) * Wrinv(i,j,k) * Xi(i,j,k) else if (grid_r(i,j,k) > (TOV_surface + rprec)) then drXi(i,j,k) = 0 drrXi(i,j,k) = 0 Xi(i,j,k) = 0 Pi(i,j,k) = 0 Xidot(i,j,k) = 0 Pidot(i,j,k) = 0 Xeta(i,j,k) = 0 end if end do end do end do !$OMP END PARALLEL DO
I attached the specific files of my thorn, the parameter file, as well as some screenshots of the evolved data.
As you can see, the boundary conditions do not affect any part along the x axis. Does any one have an idea why my loop seems to "miss" some if-cases?
Thanks a lot and best regards,
Severin
Users mailing list Users@einsteintoolkit.org http://lists.einsteintoolkit.org/mailman/listinfo/users