added natural id to ReactionStepOperatorSplit argument list

Author: ghammond86

alquimia/alquimia_fortran_interface_mod.F90

@@ -208,8 +208,8 @@ subroutine ProcessCondition(pft_engine_state,condition,props,state,aux_data,stat
 
     ! take one (or more?) reaction steps in operator split mode
     interface
-      subroutine ReactionStepOperatorSplit(pft_engine_state, delta_t, props, state, aux_data, status) bind(C)
-        use, intrinsic :: iso_c_binding, only : c_ptr, c_double
+      subroutine ReactionStepOperatorSplit(pft_engine_state, delta_t, props, state, aux_data, natural_id, status) bind(C)
+        use, intrinsic :: iso_c_binding, only : c_ptr, c_double, c_int
         use AlquimiaContainers_module, only : AlquimiaSizes,AlquimiaProblemMetaData,AlquimiaProperties,&
                  AlquimiaState,AlquimiaAuxiliaryData,AlquimiaAuxiliaryOutputData, AlquimiaEngineStatus,&
                  AlquimiaGeochemicalCondition,AlquimiaEngineFunctionality
@@ -219,6 +219,7 @@ subroutine ReactionStepOperatorSplit(pft_engine_state, delta_t, props, state, au
         type(AlquimiaProperties) :: props
         type(AlquimiaState) :: state
         type(AlquimiaAuxiliaryData) :: aux_data
+        integer(c_int),value :: natural_id
         type(AlquimiaEngineStatus) :: status
       end subroutine
     end interface
@@ -310,8 +311,8 @@ subroutine Alquimia_Fortran_ProcessCondition(this,pft_engine_state,condition,pro
       call engine_ProcessCondition(pft_engine_state,condition,props,state,aux_data,status)
     end subroutine
 
-  subroutine Alquimia_Fortran_ReactionStepOperatorSplit(this,pft_engine_state, delta_t, props, state, aux_data, status)
-    use, intrinsic :: iso_c_binding, only : c_ptr, c_double,c_f_procpointer
+  subroutine Alquimia_Fortran_ReactionStepOperatorSplit(this,pft_engine_state, delta_t, props, state, aux_data, natural_id, status)
+    use, intrinsic :: iso_c_binding, only : c_ptr, c_double,c_f_procpointer,c_int
     use AlquimiaContainers_module, only : AlquimiaSizes,AlquimiaProblemMetaData,AlquimiaProperties,&
              AlquimiaState,AlquimiaAuxiliaryData,AlquimiaAuxiliaryOutputData, AlquimiaEngineStatus,&
              AlquimiaGeochemicalCondition,AlquimiaEngineFunctionality
@@ -322,12 +323,13 @@ subroutine Alquimia_Fortran_ReactionStepOperatorSplit(this,pft_engine_state, del
     type(AlquimiaProperties) :: props
     type(AlquimiaState) :: state
     type(AlquimiaAuxiliaryData) :: aux_data
+    integer(c_int),value :: natural_id
     type(AlquimiaEngineStatus) :: status
 
     procedure(ReactionStepOperatorSplit), pointer :: engine_ReactionStepOperatorSplit
 
     call c_f_procpointer(this%c_interface%ReactionStepOperatorSplit,engine_ReactionStepOperatorSplit)
-    call engine_ReactionStepOperatorSplit(pft_engine_state, delta_t, props, state, aux_data, status)
+    call engine_ReactionStepOperatorSplit(pft_engine_state, delta_t, props, state, aux_data, natural_id, status)
   end subroutine
 
   subroutine Alquimia_Fortran_GetAuxiliaryOutput(this,pft_engine_state,props,state,aux_data,aux_out,status)

alquimia/alquimia_interface.h

@@ -93,6 +93,7 @@ extern "C" {
         AlquimiaProperties* props,
         AlquimiaState* state,
         AlquimiaAuxiliaryData* aux_data,
+        int natural_id, 
         AlquimiaEngineStatus* status);
 
     /* Access to user selected geochemical data for output, i.e. pH, 

alquimia/crunch_alquimia_interface.F90

@@ -602,10 +602,10 @@ end subroutine ProcessCondition
 
 ! **************************************************************************** !
 subroutine ReactionStepOperatorSplit(cf_engine_state, &
-     delta_t, properties, state, aux_data, status)
+     delta_t, properties, state, aux_data, natural_id, status)
 !  NOTE: Function signature is dictated by the alquimia API.
 
-  use, intrinsic :: iso_c_binding, only : c_ptr, c_double, c_f_pointer
+  use, intrinsic :: iso_c_binding, only : c_ptr, c_double, c_f_pointer, c_int
 
   use c_f_interface_module, only : f_c_string_ptr
 
@@ -630,6 +630,7 @@ subroutine ReactionStepOperatorSplit(cf_engine_state, &
   type (AlquimiaProperties), intent(in) :: properties
   type (AlquimiaState), intent(inout) :: state
   type (AlquimiaAuxiliaryData), intent(inout) :: aux_data
+  integer (c_int), value, intent(in) :: natural_id
   type (AlquimiaEngineStatus), intent(out) :: status
 
   ! local variables

alquimia/crunch_alquimia_interface.h

@@ -67,6 +67,7 @@ extern "C" {
       AlquimiaProperties* properties,
       AlquimiaState* state,
       AlquimiaAuxiliaryData* aux_data,
+      int natural_id,
       AlquimiaEngineStatus* status);
   void crunch_alquimia_getauxiliaryoutput(
       void* cf_engine_state,

alquimia/crunch_alquimia_wrappers.F90

@@ -134,6 +134,7 @@ subroutine Crunch_Alquimia_ReactionStepOperatorSplit( &
      properties, &
      state, &
      aux_data, &
+     natural_id, &
      status) bind(C)
 
   use, intrinsic :: iso_c_binding
@@ -149,10 +150,11 @@ subroutine Crunch_Alquimia_ReactionStepOperatorSplit( &
   type (AlquimiaProperties), intent(in) :: properties
   type (AlquimiaState), intent(inout) :: state
   type (AlquimiaAuxiliaryData), intent(inout) :: aux_data
+  integer (c_int), value, intent(in) :: natural_id
   type (AlquimiaEngineStatus), intent(out) :: status
 
   call ReactionStepOperatorSplit(cf_engine_state, delta_t, &
-       properties, state, aux_data, status)
+       properties, state, aux_data, natural_id, status)
 
 end subroutine Crunch_Alquimia_ReactionStepOperatorSplit
 

alquimia/pflotran_alquimia_interface.F90

@@ -472,12 +472,12 @@ end subroutine ProcessCondition
 
 ! **************************************************************************** !
 subroutine ReactionStepOperatorSplit(pft_engine_state, &
-     delta_t, properties, state, aux_data, status)
+     delta_t, properties, state, aux_data, natural_id, status)
 !  NOTE: Function signature is dictated by the alquimia API.
 
 #include "petsc/finclude/petscsys.h"
   use petscsys
-  use, intrinsic :: iso_c_binding, only : c_ptr, c_double, c_f_pointer
+  use, intrinsic :: iso_c_binding, only : c_ptr, c_double, c_f_pointer, c_int
 
   use c_f_interface_module, only : f_c_string_ptr
 
@@ -494,14 +494,14 @@ subroutine ReactionStepOperatorSplit(pft_engine_state, &
   type (AlquimiaProperties), intent(in) :: properties
   type (AlquimiaState), intent(inout) :: state
   type (AlquimiaAuxiliaryData), intent(inout) :: aux_data
+  integer (c_int) :: natural_id
   type (AlquimiaEngineStatus), intent(out) :: status
 
   ! local variables
   type(PFLOTRANEngineState), pointer :: engine_state
   PetscReal :: porosity, volume, vol_frac_prim
   PetscReal, allocatable :: guess(:)
   PetscInt :: i, num_newton_iterations, ierror
-  PetscInt, parameter :: natural_id = -999
   PetscInt, parameter :: phase_index = 1
   logical, parameter :: copy_auxdata = .true.
   class(reaction_rt_type), pointer :: reaction

alquimia/pflotran_alquimia_interface.h

@@ -67,6 +67,7 @@ extern "C" {
       AlquimiaProperties* material_properties,
       AlquimiaState* state,
       AlquimiaAuxiliaryData* aux_data,
+      int natural_id,
       AlquimiaEngineStatus* status);
   void pflotran_alquimia_getauxiliaryoutput(
       void* pft_engine_state,

alquimia/pflotran_alquimia_wrappers.F90

@@ -133,6 +133,7 @@ subroutine PFloTran_Alquimia_ReactionStepOperatorSplit( &
      properties, &
      state, &
      aux_data, &
+     natural_id, &
      status) bind(C)
 
   use, intrinsic :: iso_c_binding
@@ -148,10 +149,11 @@ subroutine PFloTran_Alquimia_ReactionStepOperatorSplit( &
   type (AlquimiaProperties), intent(in) :: properties
   type (AlquimiaState), intent(inout) :: state
   type (AlquimiaAuxiliaryData), intent(inout) :: aux_data
+  integer(c_int), value, intent(in) :: natural_id
   type (AlquimiaEngineStatus), intent(out) :: status
 
   call ReactionStepOperatorSplit(pft_engine_state, delta_t, &
-       properties, state, aux_data, status)
+       properties, state, aux_data, natural_id, status)
 
 end subroutine PFloTran_Alquimia_ReactionStepOperatorSplit
 

drivers/BatchChemDriver.c

@@ -555,6 +555,7 @@ int BatchChemDriver_Run(BatchChemDriver* driver)
     printf("BatchChemDriver: running %s\n", driver->description);
 
   double dt = driver->dt;
+  int natural_id;
 
   // Initialize the chemistry state and set up the solution vector.
   int status = BatchChemDriver_Initialize(driver);
@@ -569,10 +570,12 @@ int BatchChemDriver_Run(BatchChemDriver* driver)
       printf("BatchChemDriver: step %d (t = %g, dt = %g)\n", driver->step, driver->time, dt);
 
     // Do the chemistry step.
+    natural_id = -999;
     driver->chem.ReactionStepOperatorSplit(&driver->chem_engine,
                                            dt, &driver->chem_properties,
                                            &driver->chem_state,
                                            &driver->chem_aux_data,
+                                           natural_id,
                                            &driver->chem_status);
     if (driver->chem_status.error != 0)
     {

drivers/TransportDriver.c

@@ -669,6 +669,7 @@ static int Run_OperatorSplit(TransportDriver* driver)
   double dx = (driver->x_max - driver->x_min) / driver->num_cells;
   double dt = Min(driver->dt, driver->cfl * dx / driver->vx);
   int num_primary = driver->chem_sizes.num_primary;
+  int natural_id;
 
   // Initialize the chemistry state in each cell, and set up the solution vector.
   int status = TransportDriver_Initialize(driver);
@@ -710,10 +711,12 @@ static int Run_OperatorSplit(TransportDriver* driver)
       }
 
       // Do the chemistry step, using the advected state as input.
+      natural_id = -999;
       driver->chem.ReactionStepOperatorSplit(&driver->chem_engine,
                                              dt, &driver->chem_properties[i],
                                              &driver->advected_chem_state,
                                              &driver->advected_chem_aux_data,
+                                             natural_id,
                                              &driver->chem_status);
       if (driver->chem_status.error != 0)
       {