Merge pull request #88 from smolins/new-gases

Upgrade to interface to handle gases explicitly
LBL-EESA · Sep 11, 2023 · 3be3d1b · 3be3d1b
2 parents 043fe0e + 157a986
commit 3be3d1b
Show file tree

Hide file tree

Showing 17 changed files with 464 additions and 15 deletions.
diff --git a/alquimia/alquimia_constants.c b/alquimia/alquimia_constants.c
@@ -45,6 +45,7 @@ const char* kAlquimiaStringPH = "pH";
 const char* kAlquimiaStringMineral = "mineral";
 const char* kAlquimiaStringGas = "gas";
 const char* kAlquimiaStringCharge = "charge";
+const char* kAlquimiaStringTotalGaseous = "total_gaseous";
 
 /* Error Codes */
 const int kAlquimiaNoError = 0;

diff --git a/alquimia/alquimia_containers.F90 b/alquimia/alquimia_containers.F90
@@ -100,6 +100,8 @@ module AlquimiaContainers_module
      integer (c_int) :: num_surface_sites
      integer (c_int) :: num_ion_exchange_sites
      integer (c_int) :: num_isotherm_species
+     integer (c_int) :: num_gases
+!!     integer (c_int) :: num_gas_species
      integer (c_int) :: num_aux_integers
      integer (c_int) :: num_aux_doubles
   end type AlquimiaSizes
@@ -115,6 +117,8 @@ module AlquimiaContainers_module
      type (AlquimiaVectorDouble) :: mineral_specific_surface_area
      type (AlquimiaVectorDouble) :: surface_site_density
      type (AlquimiaVectorDouble) :: cation_exchange_capacity
+!!     type (AlquimiaVectorDouble) :: total_gas
+     type (AlquimiaVectorDouble) :: gas_concentration
   end type AlquimiaState
 
   type, public, bind(c) :: AlquimiaProperties
@@ -160,6 +164,7 @@ module AlquimiaContainers_module
      type (AlquimiaVectorString) :: ion_exchange_names
      type (AlquimiaVectorString) :: isotherm_species_names
      type (AlquimiaVectorString) :: aqueous_kinetic_names
+     type (AlquimiaVectorString) :: gas_names
   end type AlquimiaProblemMetaData
 
   type, public, bind(c) :: AlquimiaAuxiliaryOutputData
@@ -171,6 +176,7 @@ module AlquimiaContainers_module
      type (AlquimiaVectorDouble) :: primary_activity_coeff
      type (AlquimiaVectorDouble) :: secondary_free_ion_concentration
      type (AlquimiaVectorDouble) :: secondary_activity_coeff
+     type (AlquimiaVectorDouble) :: gas_partial_pressure
   end type AlquimiaAuxiliaryOutputData
 
   type, public, bind(c) :: AlquimiaAqueousConstraint

diff --git a/alquimia/alquimia_containers.h b/alquimia/alquimia_containers.h
@@ -70,6 +70,8 @@ extern "C" {
     int num_surface_sites;
     int num_ion_exchange_sites;
     int num_isotherm_species;
+    //    int num_total_gases;
+    int num_gases;
     int num_aux_integers;
     int num_aux_doubles;
   } AlquimiaSizes;
@@ -84,7 +86,9 @@ extern "C" {
     AlquimiaVectorDouble mineral_volume_fraction;  /* [-] */
     AlquimiaVectorDouble mineral_specific_surface_area; /* [m^2 mineral/m^3 bulk] */
     AlquimiaVectorDouble surface_site_density;  /* [moles/m^3 bulk] */
-    AlquimiaVectorDouble cation_exchange_capacity;  /* [moles/m^3 bulk] */
+    AlquimiaVectorDouble cation_exchange_capacity; /* [moles/m^3 bulk] */
+    //    AlquimiaVectorDouble total_gas;  /* [moles/L gas] */
+    AlquimiaVectorDouble gas_concentration;  /* [moles/L gas] */
   } AlquimiaState;
 
   typedef struct {
@@ -129,6 +133,7 @@ extern "C" {
     AlquimiaVectorString ion_exchange_names;
     AlquimiaVectorString isotherm_species_names;
     AlquimiaVectorString aqueous_kinetic_names;
+    AlquimiaVectorString gas_names;
   } AlquimiaProblemMetaData;
 
   typedef struct {
@@ -140,6 +145,7 @@ extern "C" {
     AlquimiaVectorDouble primary_activity_coeff; /* [-] */
     AlquimiaVectorDouble secondary_free_ion_concentration; /* [molality] */
     AlquimiaVectorDouble secondary_activity_coeff; /* [-] */
+    AlquimiaVectorDouble gas_partial_pressure; /* [-] */
   } AlquimiaAuxiliaryOutputData;
 
   /* 

diff --git a/alquimia/alquimia_memory.c b/alquimia/alquimia_memory.c
@@ -164,6 +164,12 @@ void AllocateAlquimiaState(const AlquimiaSizes* const sizes,
 
   AllocateAlquimiaVectorDouble(sizes->num_minerals,
                                &(state->mineral_specific_surface_area));
+
+  //  AllocateAlquimiaVectorDouble(sizes->num_gases,
+  //                             &(state->total_gas));
+
+  AllocateAlquimiaVectorDouble(sizes->num_gases,
+                               &(state->gas_concentration));
 }  /* end AllocateAlquimiaState() */
 
 void FreeAlquimiaState(AlquimiaState* state) {
@@ -174,6 +180,7 @@ void FreeAlquimiaState(AlquimiaState* state) {
     FreeAlquimiaVectorDouble(&(state->mineral_specific_surface_area));
     FreeAlquimiaVectorDouble(&(state->cation_exchange_capacity));
     FreeAlquimiaVectorDouble(&(state->surface_site_density));
+    FreeAlquimiaVectorDouble(&(state->gas_concentration));
   }
 }  /* end FreeAlquimiaState() */
 
@@ -261,6 +268,9 @@ void AllocateAlquimiaProblemMetaData(const AlquimiaSizes* const sizes,
   AllocateAlquimiaVectorString(sizes->num_aqueous_kinetics,
                                &(meta_data->aqueous_kinetic_names));
 
+  AllocateAlquimiaVectorString(sizes->num_gases,
+                               &(meta_data->gas_names));
+
 }  /* end AllocateAlquimiaProblemMetaData() */
 
 void FreeAlquimiaProblemMetaData(AlquimiaProblemMetaData* meta_data) {
@@ -273,6 +283,7 @@ void FreeAlquimiaProblemMetaData(AlquimiaProblemMetaData* meta_data) {
     FreeAlquimiaVectorString(&(meta_data->ion_exchange_names));
     FreeAlquimiaVectorString(&(meta_data->isotherm_species_names));
     FreeAlquimiaVectorString(&(meta_data->aqueous_kinetic_names));
+    FreeAlquimiaVectorString(&(meta_data->gas_names));
   }
 }  /* end FreeAlquimiaProblemMetaData() */
 
@@ -303,6 +314,9 @@ void AllocateAlquimiaAuxiliaryOutputData(const AlquimiaSizes* const sizes,
                                &(aux_output->secondary_free_ion_concentration));
   AllocateAlquimiaVectorDouble(sizes->num_aqueous_complexes,
                                &(aux_output->secondary_activity_coeff));
+
+  AllocateAlquimiaVectorDouble(sizes->num_gases,
+                               &(aux_output->gas_partial_pressure));
 
 }  /* end AllocateAlquimiaAuxiliaryOutputData() */
 
@@ -315,6 +329,7 @@ void FreeAlquimiaAuxiliaryOutputData(AlquimiaAuxiliaryOutputData* aux_output) {
     FreeAlquimiaVectorDouble(&(aux_output->primary_activity_coeff));
     FreeAlquimiaVectorDouble(&(aux_output->secondary_free_ion_concentration));
     FreeAlquimiaVectorDouble(&(aux_output->secondary_activity_coeff));
+	FreeAlquimiaVectorDouble(&(aux_output->gas_partial_pressure));
   }
 }  /* end FreeAlquimiaAuxiliaryOutputData() */
 

diff --git a/alquimia/alquimia_util.c b/alquimia/alquimia_util.c
@@ -131,6 +131,7 @@ void CopyAlquimiaProblemMetaData(const AlquimiaProblemMetaData* const source,
   CopyAlquimiaVectorString(&source->ion_exchange_names, &destination->ion_exchange_names);
   CopyAlquimiaVectorString(&source->isotherm_species_names, &destination->isotherm_species_names);
   CopyAlquimiaVectorString(&source->aqueous_kinetic_names, &destination->aqueous_kinetic_names);
+  CopyAlquimiaVectorString(&source->gas_names, &destination->gas_names);
 }
 
 void CopyAlquimiaProperties(const AlquimiaProperties* const source, 
@@ -164,6 +165,7 @@ void CopyAlquimiaState(const AlquimiaState* const source,
   CopyAlquimiaVectorDouble(&source->mineral_specific_surface_area, &destination->mineral_specific_surface_area);
   CopyAlquimiaVectorDouble(&source->surface_site_density, &destination->surface_site_density);
   CopyAlquimiaVectorDouble(&source->cation_exchange_capacity, &destination->cation_exchange_capacity);
+  CopyAlquimiaVectorDouble(&source->gas_concentration, &destination->gas_concentration);
 }
 
 void CopyAlquimiaAuxiliaryData(const AlquimiaAuxiliaryData* const source, 
@@ -184,6 +186,7 @@ void CopyAlquimiaAuxiliaryOutputData(const AlquimiaAuxiliaryOutputData* const so
   CopyAlquimiaVectorDouble(&source->primary_activity_coeff, &destination->primary_activity_coeff);
   CopyAlquimiaVectorDouble(&source->secondary_free_ion_concentration, &destination->secondary_free_ion_concentration);
   CopyAlquimiaVectorDouble(&source->secondary_activity_coeff, &destination->secondary_activity_coeff);
+  CopyAlquimiaVectorDouble(&source->gas_partial_pressure, &destination->gas_partial_pressure);
 }
 
 void CopyAlquimiaGeochemicalCondition(const AlquimiaGeochemicalCondition* const source, 
@@ -397,6 +400,7 @@ void PrintAlquimiaSizes(const AlquimiaSizes* const sizes, FILE* file) {
   fprintf(file, "     num aqueous kinetics : %d\n", sizes->num_aqueous_kinetics);
   fprintf(file, "     num surface sites : %d\n", sizes->num_surface_sites);
   fprintf(file, "     num ion exchange sites : %d\n", sizes->num_ion_exchange_sites);
+  fprintf(file, "     num gases: %d\n", sizes->num_gases);
   fprintf(file, "     num auxiliary integers : %d\n", sizes->num_aux_integers);
   fprintf(file, "     num auxiliary doubles : %d\n", sizes->num_aux_doubles);
 }  /* end PrintAlquimiaSizes() */
@@ -423,6 +427,7 @@ void PrintAlquimiaProblemMetaData(const AlquimiaProblemMetaData* const meta_data
   PrintAlquimiaVectorString("ion exchange names", &(meta_data->ion_exchange_names), file);
   PrintAlquimiaVectorString("isotherm species names", &(meta_data->isotherm_species_names), file);
   PrintAlquimiaVectorString("aqueous kinetic names", &(meta_data->aqueous_kinetic_names), file);
+  PrintAlquimiaVectorString("gas names", &(meta_data->gas_names), file);
 }  /* end PrintAlquimiaProblemMetaData() */
 
 void PrintAlquimiaProperties(const AlquimiaProperties* const mat_prop, FILE* file) {
@@ -455,6 +460,7 @@ void PrintAlquimiaState(const AlquimiaState* const state, FILE* file) {
                             &(state->cation_exchange_capacity), file);
   PrintAlquimiaVectorDouble("surface_site_density",
                             &(state->surface_site_density), file);
+  PrintAlquimiaVectorDouble("gas_concentration", &(state->gas_concentration), file);
 }  /* end PrintAlquimiaState() */
 
 void PrintAlquimiaAuxiliaryData(const AlquimiaAuxiliaryData* const aux_data, FILE* file) {
@@ -481,6 +487,8 @@ void PrintAlquimiaAuxiliaryOutputData(const AlquimiaAuxiliaryOutputData* const a
                             &(aux_output->secondary_free_ion_concentration), file);
   PrintAlquimiaVectorDouble("secondary activity coeff",
                             &(aux_output->secondary_activity_coeff), file);
+  PrintAlquimiaVectorDouble("gas_partial_pressure",
+                            &(aux_output->gas_partial_pressure), file);
 }  /* end PrintAlquimiaAuxiliaryOutputData() */
 
 void PrintAlquimiaGeochemicalConditionVector(const AlquimiaGeochemicalConditionVector* const condition_list, FILE* file) {