Added support for porosity and forcing functions to registry.

TDycores-Project · Aug 5, 2021 · 130a002 · 130a002
1 parent 6e702a0
commit 130a002
Show file tree

Hide file tree

Showing 2 changed files with 91 additions and 30 deletions.
diff --git a/demo/transient/transient_snes_mpfaof90.F90 b/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
 
@@ -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
 
@@ -163,6 +163,8 @@ program main
 
   ! 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);
@@ -310,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
@@ -333,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);

diff --git a/src/f90-mod/tdycoremod.F90 b/src/f90-mod/tdycoremod.F90
@@ -448,6 +448,17 @@ subroutine TDyFunction(tdy, x, f, dummy, 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)
@@ -481,7 +492,7 @@ function RegisterFn(name, func) bind (c, name="TDyRegisterFunction") result(ierr
     ierr = RegisterFn(FtoCString(name), c_funloc(func))
   end subroutine
 
-  subroutine TDySelectBoundaryPressureFn(tdy, name, ierr)
+  subroutine TDySelectPorosityFunction(tdy, name, ierr)
     use, intrinsic :: iso_c_binding
     use tdycoredef
     implicit none
@@ -492,17 +503,73 @@ subroutine TDySelectBoundaryPressureFn(tdy, name, ierr)
     type(c_funptr)                  :: c_func
     procedure(TDyFunction), pointer :: f_func
 
-    interface
-      function GetFn(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
+    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 = GetFn(FtoCString(name), c_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
@@ -518,17 +585,7 @@ subroutine TDySelectBoundaryVelocityFn(tdy, name, ierr)
     type(c_funptr)                  :: c_func
     procedure(TDyFunction), pointer :: f_func
 
-    interface
-      function GetFn(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
-
-    ierr = GetFn(FtoCString(name), c_func)
+    ierr = GetRegFn(FtoCString(name), c_func)
     call c_f_procpointer(c_func, f_func)
     call TDySetBoundaryVelocityFn(tdy, f_func, 0, ierr)
   end subroutine