nferraro/fusion-io

Libraries providing common interface to data from various fusion codes

fusion-io's Issues

Add license

Preferably MIT.

Changes to fpy.py

@andreas-kleiner I'm thinking of making a few changes to how fpy.py works. Let me know your thoughts.

1. Store the absolute path as sim_data.filename, that way it can be accessed if the filename is required. This would mean that we don't have to pass around sim_data and filename separately and it would ensure consistency.

2. Initialize sim_data.hint to None, then instead of calling sim_data.hint directly, we do sim_data.get_hint(). That function will return sim_data.hint, but if it's None, it will call fio_py.allocate_hint(). By moving this outside of __init__ we could initialize the sim_data object and use it for anything that doesn't require a hint without waiting for the mesh connectivity to be computed. That would only get done when it's used.

Latest Fortran interface crashes in fio_eval_field_f

I encountered a problem with the latest version of the Fortran interface on the Greene cluster. When trying to read a field with fit_eval_field_f, the my code crashes.

I think the problem is that the C-version of fio_eval_field is called with a constant (0):

fusion-io/fusion_io/fusion_io_fortran.F90

Lines 40 to 54 in 4ba28ce

	subroutine fio_eval_field_f(ifield, x, v, ierr)
	integer, intent(in) :: ifield
	real, intent(in), dimension(*) :: x
	real, intent(out), dimension(*) :: v
	integer, intent(out) :: ierr
	call fio_eval_field(ifield, x, v, 0, ierr)
	end subroutine fio_eval_field_f
	subroutine fio_eval_field_deriv_f(ifield, x, v, ierr)
	integer, intent(in) :: ifield
	real, intent(in), dimension(*) :: x
	real, intent(out), dimension(*) :: v
	integer, intent(out) :: ierr
	call fio_eval_field_deriv(ifield, x, v, 0, ierr)
	end subroutine fio_eval_field_deriv_f

I modified the code slightly by declaring an integer variable s and setting s=0

  subroutine fio_eval_field_f(ifield, x, v, ierr)
    integer, intent(in) :: ifield
    real, intent(in), dimension(*) :: x
    real, intent(out), dimension(*) :: v
    integer, intent(out) :: ierr
    integer :: s
    s=0
    call fio_eval_field(ifield, x, v, s, ierr)
  end subroutine fio_eval_field_f

  subroutine fio_eval_field_deriv_f(ifield, x, v, ierr)
    integer, intent(in) :: ifield
    real, intent(in), dimension(*) :: x
    real, intent(out), dimension(*) :: v
    integer, intent(out) :: ierr
    integer :: s
    s=0
    call fio_eval_field_deriv(ifield, x, v, s, ierr)
  end subroutine fio_eval_field_deriv_f

This turned out to work. The only way I can make sense of this is that somewhere down the call-tree, a program tries to write into s, which results in a segfault if s is a constant.

Return 0 for points out-of-bounds

Currently, the python interface crashes if you ask for the value of a field at a point outside of the domain. As you don't always know the domain a priori, or you want to look at things on a rectangular grid that the domain boundary may run through, I would suggest that zero gets returned instead of null. @nferraro If you could give me the ability to write to this repository (at least new branches if not master), I can push the suggested changes for your review. Alternatively, I've included the suggested changes below that worked for me.

diff --git a/fusion_io/python_interface.cpp b/fusion_io/python_interface.cpp
index 338170e..e3e4c41 100644
--- a/fusion_io/python_interface.cpp
+++ b/fusion_io/python_interface.cpp
@@ -148,6 +148,8 @@ PyObject* fio_eval_scalar_field_py(PyObject* self, PyObject *args)
 
   double v;
   int ierr = fio_eval_field(ifield, x, &v);
+  if(ierr == FIO_OUT_OF_BOUNDS)
+     return Py_BuildValue("d", 0.);
   if(ierr != FIO_SUCCESS)
     return NULL; 
 
@@ -180,6 +182,8 @@ PyObject* fio_eval_vector_field_py(PyObject* self, PyObject *args)
 
   double v[3];
   int ierr = fio_eval_field(ifield, x, v);
+  if(ierr == FIO_OUT_OF_BOUNDS)
+     return Py_BuildValue("(ddd)", 0., 0., 0.);
   if(ierr != FIO_SUCCESS)
     return NULL;