resetting exchanger and surfaces from solutions in phreeqcRM

[oatteia]

Hi,
I work on a coupled solution using OpenFoam for transport and phreeqcRM for chemistry. Unfortunately when you work on a large problem, the model can crash in the middle and I want to be able to restart the model. When restart, I am able to retrieve the composition of the solution in any cells, and I would like to start a new simulation by equilibrating my exchanger and surfaces with these solutions. I know that it is possible to initialize phreeqc with solutions but it would mean building 100000 solutions to equilibrate. Is there a more direct way : provide the cell chemistry and ask phreeqcRM to equilibrate exchanger and surfaces with these dissolved concentrations cells (I did not find the function)
thanks again for phreeqc!
Olivier

[oatteia]

Me again,
I found the beginning of a solution: it is possible to dump all phreeqcRM info to a file, in raw format, which works...but I cannot find the way to retrieve it in a new simulation (when phreeqc fails, it breaks it all, so i cannot make refernce to a state)
olivier

[dlparkhurst]

Yes, DumpModule is intended to be a way to restart a simulation.

If a simulation fails and the program aborts, then any SaveState information will be lost. So, you can either capture the situation with a C++ try/catch construct and try to recover from a saved state, or you can write a way to restart from the last use of DumpModule.

PHAST has a restart mechanism that you could look at and decide whether you want to implement it. Here is the code that writes the restart file:

Code: [Select]

/* ---------------------------------------------------------------------- */
IRM_RESULT
FileHandler::WriteRestart(int *print_restart)
/* ---------------------------------------------------------------------- */
{
PhreeqcRM * Reaction_module_ptr = PhreeqcRM::GetInstance(this->rm_id);
if (Reaction_module_ptr)
{
//Reaction_module_ptr->GetSaturation(this->saturation);
int mpi_myself = Reaction_module_ptr->GetMpiMyself();
if (print_restart != 0)
{
gzFile restart_file;
#ifdef USE_GZ
std::string temp_name("temp_restart_file.gz");
#else
std::string temp_name("temp_restart_file");
#endif

std::string name(Reaction_module_ptr->GetFilePrefix());
std::string backup_name(name);
if (mpi_myself == 0)
{
#ifdef USE_GZ
name.append(".restart.gz");
backup_name.append(".restart.gz~");
#else
name.append(".restart");
backup_name.append(".restart~");
#endif

// open file
restart_file = gzopen(temp_name.c_str(), "wb");
if (restart_file == NULL)
{
std::ostringstream errstr;
errstr << "Temporary restart file could not be opened: " << temp_name;
error_msg(errstr.str().c_str(), 1);
throw PhreeqcRMStop();
}

// write header
int count_chemistry = Reaction_module_ptr->GetChemistryCellCount();
gzprintf(restart_file, "%s\n", "#PHAST restart file");
gzprintf(restart_file, "%s%s\n", "#Prefix: ", Reaction_module_ptr->GetFilePrefix().c_str());
time_t now = ::time(NULL);
gzprintf(restart_file, "%s%s\n", "#Date: ", ctime(&now));
gzprintf(restart_file, "%s%e\n", "#Current model time: ", Reaction_module_ptr->GetTime());
gzprintf(restart_file, "%s%d\n", "#Grid cells: ", Reaction_module_ptr->GetGridCellCount());
gzprintf(restart_file, "%s%d\n", "#Chemistry cells: ", count_chemistry);

// write index
gzprintf(restart_file, "%d\n", count_chemistry );
for (int j = 0; j < count_chemistry; j++) /* j is count_chem number */
{
int i = Reaction_module_ptr->GetBackwardMapping()[j][0]; /* i is nxyz number */
gzprintf(restart_file, "%e %e %e %d ",  x_node[i], y_node[i], z_node[i], j);
// write all solutions
//gzprintf(restart_file, "%d ", this->ic[7 * i]);
// solution, use -1 if cell is dry
//if (this->saturation[i] > 0.0)
//{
// gzprintf(restart_file, "%d ", this->ic[7 * i]);
//}
//else
//{
// gzprintf(restart_file, "-1 ");
//}
if (this->saturation[i] >= 1.0)
{
gzprintf(restart_file, "%d ", 1);
}
else if (this->saturation[i] > 0.0)
{
gzprintf(restart_file, "%d ", 0);
}
else
{
gzprintf(restart_file, "-1 ");
}
// pp_assemblage
gzprintf(restart_file, "%d ", this->ic[7 * i + 1]);
// exchange
gzprintf(restart_file, "%d ", this->ic[7 * i + 2]);
// surface
gzprintf(restart_file, "%d ", this->ic[7 * i + 3]);
// gas_phase
gzprintf(restart_file, "%d ", this->ic[7 * i + 4]);
// solid solution
gzprintf(restart_file, "%d ", this->ic[7 * i + 5]);
// kinetics
gzprintf(restart_file, "%d \n", this->ic[7 * i + 6]);
}
gzclose(restart_file);
}
// write data
Reaction_module_ptr->SetDumpFileName(temp_name.c_str());
Reaction_module_ptr->DumpModule(true, true);

// rename files
if (mpi_myself == 0)
{
PhreeqcRM::FileRename(temp_name.c_str(), name.c_str(), backup_name.c_str());
}
return IRM_OK;
}
return IRM_INVALIDARG;
}
return IRM_BADINSTANCE;
}

The code to process restart files is more complicated. There is an input processor that converts flow and transport files into a temporary file that is actually read by the model. I'll just say that the initial conditions can be set by processing the restart file. The actual input for the model looks something like this:

Code: [Select]

Line 0: CHEMISTRY_IC
Line 1: -box 0.0 0.0 0.0 1000.0 100.0 10.0
Line 8: -solid_solutions restart project.restart.gz

Input for PHAST is grid independent, so a volume is defined for which the restart file provides the spatial distribution, in this case for the solid solution initial conditions. The restart file, which may be from a different grid than the current model, is read with X, Y, Z, and solid solution composition. The solid solution compositions are interpolated to the current grid; I think the solid solution composition of the point in the restart file that is closest to the grid point is used. If you want one restart file to populate the entire model, then the volume is simply the domain and every reactant -solution, -equilibrium_phases, etc use the same restart file.

[oatteia]

Hi me again on the same subject.
looked at the dump module, I already have the way to store all data. But I did not find any class in phreeqcRM to repopulate all cells minerals surfaces or exchanger (which are different in any cell). Can you detail how PHAST can change the values of minerals concentrations (for surface and exchanger it might be possible to iteratively equilibrate the restart solutions with the surface and exchangers)
sicnerely
Olivier

[dlparkhurst]

I'm not sure I entirely understand your question, but I think it is how to save a state and restart the calculation.

If you have dumped the contents of all the cells to a file in the raw formats (including solution, equilibrium phases, surface, exchange, and others), you would need to read the entire file into the initial phreeqc instance, and then distribute the solutions, equilibrium phases, surfaces, and exchangers (and any other reactants) to the workers with InitialPhreeqc2Module method.

You will also have to consider how to save any boundary condition data that are in the initial phreeqc instance at the point of dumping the cells. I think you may need to dump the initial phreeqc instance solutions that are used as boundary conditions to a separate file using Phreeqc's DUMP keyword. Once you have reset the cells in the workers, you could read the dump file from the initial phreeqc instance. Alternatively, you could reuse the file that was used to initialize the calculation that has the boundary condition definitions.

The approach is a resource hog because, at one point, you will have all of the data in the initial phreeqc instance and in the workers, so double the memory. As soon as you have transferred to the workers, you can cleanup the initial phreeqc instance (DELETE keyword) and read the boundary conditions back into the initial phreeqc instance.