Merge pull request #192 from TDycores-Project/jeff-cohere/named_bcs

Added a function registry for referring to functions by name.

demo/transient/transient_snes_mpfaof90.F90

@@ -9,8 +9,8 @@ module snes_mpfaof90mod
    subroutine PorosityFunction(tdy,x,theta,dummy,ierr)
      implicit none
      TDy                    :: tdy
-     PetscReal, intent(in) :: x
-     PetscReal, intent(out):: theta
+     PetscReal, intent(in)  :: x(3)
+     PetscReal, intent(out) :: theta
      integer                :: dummy(*)
      PetscErrorCode :: ierr
 
@@ -49,8 +49,8 @@ end subroutine ResidualSaturation
   subroutine PorosityFunctionPFLOTRAN(tdy,x,theta,dummy,ierr)
     implicit none
     TDy                    :: tdy
-    PetscReal, intent(in) :: x
-    PetscReal, intent(out):: theta
+    PetscReal, intent(in)  :: x(3)
+    PetscReal, intent(out) :: theta
     integer                :: dummy(*)
     PetscErrorCode :: ierr
 
@@ -84,13 +84,13 @@ subroutine MaterialPropAlpha_PFLOTRAN(alpha)
     PetscReal, intent(out) :: alpha
     alpha = 1.d-4
   end subroutine MaterialPropAlpha_PFLOTRAN
-    
+
   subroutine MaterialPropM_PFLOTRAN(m)
     implicit none
     PetscReal, intent(out) :: m
     m = 0.3d0
   end subroutine MaterialPropM_PFLOTRAN
-    
+
   subroutine ResidualSat_PFLOTRAN(resSat)
     implicit none
     PetscReal resSat
@@ -100,8 +100,8 @@ end subroutine ResidualSat_PFLOTRAN
   subroutine PressureFunction(tdy,x,pressure,dummy,ierr)
     implicit none
     TDy                    :: tdy
-    PetscReal, intent(in) :: x(3)
-    PetscReal, intent(out):: pressure
+    PetscReal, intent(in)  :: x(3)
+    PetscReal, intent(out) :: pressure
     integer                :: dummy(*)
     PetscErrorCode :: ierr
 
@@ -160,6 +160,13 @@ program main
 
   call TDyInit(ierr);
   CHKERRA(ierr);
+
+  ! Register some functions.
+  call TDyRegisterFunction("p0", PressureFunction, ierr)
+  call TDyRegisterFunction("porosity", PorosityFunction, ierr)
+  call TDyRegisterFunction("porosity_pflotran", PorosityFunctionPFLOTRAN, ierr)
+  CHKERRA(ierr);
+
   call TDyCreate(tdy, ierr);
   CHKERRA(ierr);
   call TDySetDiscretizationMethod(tdy,MPFA_O,ierr);
@@ -305,7 +312,8 @@ program main
   CHKERRA(ierr);
 
   if (pflotran_consistent) then
-     call TDySetPorosityFunction(tdy,PorosityFunctionPFLOTRAN,0,ierr);
+!     call TDySetPorosityFunction(tdy,PorosityFunctionPFLOTRAN,0,ierr);
+     call TDySelectPorosityFunction(tdy,"porosity_pflotran",ierr);
      CHKERRA(ierr);
 
      do c = 1,ncell
@@ -328,7 +336,8 @@ program main
 
   else
 
-     call TDySetPorosityFunction(tdy,PorosityFunction,0,ierr);
+!     call TDySetPorosityFunction(tdy,PorosityFunction,0,ierr);
+     call TDySelectPorosityFunction(tdy,"porosity",ierr);
      CHKERRA(ierr);
 
      call TDySetBlockPermeabilityValuesLocal(tdy,ncell,index,blockPerm,ierr);
@@ -340,7 +349,8 @@ program main
   end if
 
   if (bc_type == MPFAO_DIRICHLET_BC .OR. bc_type == MPFAO_SEEPAGE_BC ) then
-     call TDySetBoundaryPressureFn(tdy,PressureFunction,0,ierr);
+!     call TDySetBoundaryPressureFn(tdy,PressureFunction,0,ierr);
+     call TDySelectBoundaryPressureFn(tdy,"p0",ierr);
      CHKERRQ(ierr)
   endif
 

include/tdycore.h

@@ -74,8 +74,9 @@ PETSC_EXTERN PetscClassId TDY_CLASSID;
 PETSC_EXTERN PetscLogEvent TDy_ComputeSystem;
 
 /* ---------------------------------------------------------------- */
-PETSC_EXTERN PetscErrorCode TDyInit(int argc, char* argv[]);
+PETSC_EXTERN PetscErrorCode TDyInit(int, char*[]);
 PETSC_EXTERN PetscErrorCode TDyInitNoArguments(void);
+PETSC_EXTERN PetscErrorCode TDyOnFinalize(void (*)(void));
 PETSC_EXTERN PetscErrorCode TDyFinalize(void);
 
 PETSC_EXTERN PetscErrorCode TDyCreate(TDy*);
@@ -102,8 +103,12 @@ PETSC_EXTERN PetscErrorCode TDySetSoilDensityFunction(TDy,PetscErrorCode(*)(TDy,
 PETSC_EXTERN PetscErrorCode TDySetSoilSpecificHeatFunction(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*);
 
 // Set boundary and source-sink: via PETSc operations
+PetscErrorCode TDyRegisterFunction(const char*, PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*));
 PETSC_EXTERN PetscErrorCode TDySetForcingFunction(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*);
 PETSC_EXTERN PetscErrorCode TDySetEnergyForcingFunction(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*);
+PETSC_EXTERN PetscErrorCode TDySelectBoundaryPressureFn(TDy,const char*,void*);
+PETSC_EXTERN PetscErrorCode TDySelectBoundaryTemperatureFn(TDy,const char*,void*);
+PETSC_EXTERN PetscErrorCode TDySelectBoundaryVelocityFn(TDy,const char*,void*);
 PETSC_EXTERN PetscErrorCode TDySetBoundaryPressureFn(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*);
 PETSC_EXTERN PetscErrorCode TDySetBoundaryTemperatureFn(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*);
 PETSC_EXTERN PetscErrorCode TDySetBoundaryVelocityFn(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*);

src/f90-mod/tdycoremod.F90

@@ -437,6 +437,28 @@ subroutine TDyPostSolveSNESSolver(a,b,z)
      end subroutine TDyPostSolveSNESSolver
   end interface
 
+  abstract interface
+    subroutine TDyFunction(tdy, x, f, dummy, ierr)
+      use tdycoredef
+      TDy :: tdy
+      PetscReal, intent(in)  :: x(3)
+      PetscReal, intent(out) :: f
+      integer                :: dummy(*)
+      PetscErrorCode         :: ierr
+    end subroutine
+  end interface
+
+  ! We use GetRegFn to retrieve function pointers from the C registry.
+  interface
+    function GetRegFn(name, c_func) bind (c, name="TDyGetFunction") result(ierr)
+      use, intrinsic :: iso_c_binding
+      implicit none
+      type(c_ptr), value :: name
+      type(c_funptr) :: c_func
+      integer(c_int) :: ierr
+    end function
+  end interface
+
   contains
 
   subroutine TDyInit(ierr)
@@ -450,4 +472,145 @@ subroutine TDyInit(ierr)
      call TDyInitNoArguments(ierr)
   end subroutine TDyInit
 
+  subroutine TDyRegisterFunction(name, func, ierr)
+    use, intrinsic :: iso_c_binding
+    implicit none
+    character(len=*), intent(in)   :: name
+    procedure(TDyFunction)         :: func
+    PetscErrorCode                 :: ierr
+
+    interface
+      function RegisterFn(name, func) bind (c, name="TDyRegisterFunction") result(ierr)
+        use, intrinsic :: iso_c_binding
+        implicit none
+        type(c_ptr), value :: name
+        type(c_funptr), value :: func
+        integer(c_int) :: ierr
+      end function
+    end interface
+
+    ierr = RegisterFn(FtoCString(name), c_funloc(func))
+  end subroutine
+
+  subroutine TDySelectPorosityFunction(tdy, name, ierr)
+    use, intrinsic :: iso_c_binding
+    use tdycoredef
+    implicit none
+    TDy :: tdy
+    character(len=*), intent(in)   :: name
+    PetscErrorCode                 :: ierr
+
+    type(c_funptr)                  :: c_func
+    procedure(TDyFunction), pointer :: f_func
+
+    ierr = GetRegFn(FtoCString(name), c_func)
+    call c_f_procpointer(c_func, f_func)
+    call TDySetPorosityFunction(tdy, f_func, 0, ierr)
+  end subroutine
+
+  subroutine TDySelectForcingFunction(tdy, name, ierr)
+    use, intrinsic :: iso_c_binding
+    use tdycoredef
+    implicit none
+    TDy :: tdy
+    character(len=*), intent(in)   :: name
+    PetscErrorCode                 :: ierr
+
+    type(c_funptr)                  :: c_func
+    procedure(TDyFunction), pointer :: f_func
+
+    ierr = GetRegFn(FtoCString(name), c_func)
+    call c_f_procpointer(c_func, f_func)
+    call TDySetForcingFunction(tdy, f_func, 0, ierr)
+  end subroutine
+
+! Uncomment this when we are ready to set energy forcing fns in Fortran.
+!  subroutine TDySelectEnergyForcingFunction(tdy, name, ierr)
+!    use, intrinsic :: iso_c_binding
+!    use tdycoredef
+!    implicit none
+!    TDy :: tdy
+!    character(len=*), intent(in)   :: name
+!    PetscErrorCode                 :: ierr
+!
+!    type(c_funptr)                  :: c_func
+!    procedure(TDyFunction), pointer :: f_func
+!
+!    ierr = GetRegFn(FtoCString(name), c_func)
+!    call c_f_procpointer(c_func, f_func)
+!    call TDySetEnergyForcingFunction(tdy, f_func, 0, ierr)
+!  end subroutine
+
+! Uncomment this when we are ready to set permeability functions programmatically.
+!  subroutine TDySelectPermeabilityFunction(tdy, name, ierr)
+!    use, intrinsic :: iso_c_binding
+!    use tdycoredef
+!    implicit none
+!    TDy :: tdy
+!    character(len=*), intent(in)   :: name
+!    PetscErrorCode                 :: ierr
+!
+!    type(c_funptr)                  :: c_func
+!    procedure(TDyFunction), pointer :: f_func
+!
+!    ierr = GetRegFn(FtoCString(name), c_func)
+!    call c_f_procpointer(c_func, f_func)
+!    call TDySetPermeabilityFunction(tdy, f_func, 0, ierr)
+!  end subroutine
+
+  subroutine TDySelectBoundaryPressureFn(tdy, name, ierr)
+    use, intrinsic :: iso_c_binding
+    use tdycoredef
+    implicit none
+    TDy :: tdy
+    character(len=*), intent(in)   :: name
+    PetscErrorCode                 :: ierr
+
+    type(c_funptr)                  :: c_func
+    procedure(TDyFunction), pointer :: f_func
+
+    ierr = GetRegFn(FtoCString(name), c_func)
+    call c_f_procpointer(c_func, f_func)
+    call TDySetBoundaryPressureFn(tdy, f_func, 0, ierr)
+  end subroutine
+
+  subroutine TDySelectBoundaryVelocityFn(tdy, name, ierr)
+    use, intrinsic :: iso_c_binding
+    use tdycoredef
+    implicit none
+    TDy :: tdy
+    character(len=*), intent(in)   :: name
+    PetscErrorCode                 :: ierr
+
+    type(c_funptr)                  :: c_func
+    procedure(TDyFunction), pointer :: f_func
+
+    ierr = GetRegFn(FtoCString(name), c_func)
+    call c_f_procpointer(c_func, f_func)
+    call TDySetBoundaryVelocityFn(tdy, f_func, 0, ierr)
+  end subroutine
+
+  ! Here's a function that converts a Fortran string to a C string and
+  ! stashes it in TDycore's Fortran string registry. This allows us to
+  ! create more expressive (and standard) Fortran interfaces.
+  function FtoCString(f_string) result(c_string)
+     use, intrinsic :: iso_c_binding
+     implicit none
+     character(len=*), target :: f_string
+     character(len=:), pointer :: f_ptr
+     type(c_ptr) :: c_string
+
+     interface
+       function NewCString(f_str_ptr, f_str_len) bind (c, name="NewCString") result(c_string)
+         use, intrinsic :: iso_c_binding
+         type(c_ptr), value :: f_str_ptr
+         integer(c_int), value :: f_str_len
+         type(c_ptr) :: c_string
+       end function NewCString
+     end interface
+
+     f_ptr => f_string
+     c_string = NewCString(c_loc(f_ptr), len(f_string))
+   end function FtoCString
+
 end module tdycore

src/interface/ftn/tdycoref.c

@@ -130,6 +130,7 @@
 #define tdydestroy_                                    tdydestroy
 #endif
 
+// Legacy Fortran function pointer registry
 static struct {
   PetscFortranCallbackId porosity;
   PetscFortranCallbackId permeability;
@@ -142,7 +143,6 @@ static struct {
 #endif
 } _cb;
 
-
 #if defined(__cplusplus)
 extern "C" {
 #endif

src/tdyconditions.c

@@ -1,9 +1,56 @@
 #include <private/tdycoreimpl.h>
+#include <petsc/private/khash/khash.h>
 
 /*
   Boundary and source-sink conditions are set by PETSc operations
 */
 
+// Here's a registry of functions that can be used for boundary conditions and
+// forcing terms.
+typedef PetscErrorCode(*Function)(TDy, PetscReal*, PetscReal*, void*);
+KHASH_MAP_INIT_STR(TDY_FUNC_MAP, Function)
+static khash_t(TDY_FUNC_MAP)* funcs_ = NULL;
+
+// This function is called on finalization to destroy the function registry.
+static void DestroyFunctionRegistry() {
+  kh_destroy(TDY_FUNC_MAP, funcs_);
+}
+
+PetscErrorCode TDyRegisterFunction(const char* name, Function f) {
+  PetscFunctionBegin;
+  if (funcs_ == NULL) {
+    funcs_ = kh_init(TDY_FUNC_MAP);
+    TDyOnFinalize(DestroyFunctionRegistry);
+  }
+
+  int retval;
+  khiter_t iter = kh_put(TDY_FUNC_MAP, funcs_, name, &retval);
+  kh_val(funcs_, iter) = f;
+  PetscFunctionReturn(0);
+}
+
+PetscErrorCode TDyGetFunction(const char* name, Function* f) {
+  PetscFunctionBegin;
+  int ierr;
+
+  if (funcs_ != NULL) {
+    khiter_t iter = kh_get(TDY_FUNC_MAP, funcs_, name);
+    if (iter != kh_end(funcs_)) { // found it!
+      *f = kh_val(funcs_, iter);
+    } else {
+      ierr = -1;
+      SETERRQ(MPI_COMM_WORLD, ierr, "Function not found!");
+      return ierr;
+    }
+  } else {
+    ierr = -1;
+    SETERRQ(MPI_COMM_WORLD, ierr, "No functions have been registered!");
+    return ierr;
+  }
+  PetscFunctionReturn(0);
+}
+
+
 PetscErrorCode TDySetForcingFunction(TDy tdy, PetscErrorCode(*f)(TDy,PetscReal*,PetscReal*,void*),void *ctx) {
   PetscFunctionBegin;
   if (f) tdy->ops->computeforcing = f;
@@ -18,6 +65,33 @@ PetscErrorCode TDySetEnergyForcingFunction(TDy tdy, PetscErrorCode(*f)(TDy,Petsc
   PetscFunctionReturn(0);
 }
 
+PetscErrorCode TDySelectBoundaryPressureFn(TDy tdy, const char* name, void* ctx) {
+  PetscFunctionBegin;
+  int ierr;
+  Function f;
+  ierr = TDyGetFunction(name, &f); CHKERRQ(ierr);
+  ierr = TDySetBoundaryPressureFn(tdy, f, ctx); CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+PetscErrorCode TDySelectBoundaryTemperatureFn(TDy tdy, const char* name, void* ctx) {
+  PetscFunctionBegin;
+  int ierr;
+  Function f;
+  ierr = TDyGetFunction(name, &f); CHKERRQ(ierr);
+  ierr = TDySetBoundaryTemperatureFn(tdy, f, ctx); CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+PetscErrorCode TDySelectBoundaryVelocityFn(TDy tdy, const char* name, void* ctx) {
+  PetscFunctionBegin;
+  int ierr;
+  Function f;
+  ierr = TDyGetFunction(name, &f); CHKERRQ(ierr);
+  ierr = TDySetBoundaryVelocityFn(tdy, f, ctx); CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
 PetscErrorCode TDySetBoundaryTemperatureFn(TDy tdy, PetscErrorCode(*f)(TDy,PetscReal*,PetscReal*,void*),void *ctx) {
   PetscFunctionBegin;
   if (f) tdy->ops->compute_boundary_temperature = f;