Saved surface for numerical derivatives

This commit is contained in:
David Parkhurst 2021-07-17 08:49:17 -06:00
parent df0d68b9d5
commit 56975a7e0c
3 changed files with 36 additions and 10 deletions

View File

@ -5560,6 +5560,7 @@ numerical_jacobian(void)
std::vector<class phase> base_phases;
double base_mass_water_bulk_x = 0, base_moles_h2o = 0;
cxxGasPhase base_gas_phase;
cxxSurface base_surface;
if (!
(numerical_deriv ||
@ -5571,6 +5572,10 @@ numerical_jacobian(void)
calculating_deriv = TRUE;
//jacobian_sums();
if (use.Get_surface_ptr() != NULL)
{
base_surface = *use.Get_surface_ptr();
}
if (use.Get_gas_phase_ptr() != NULL)
{
//cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
@ -5761,6 +5766,10 @@ numerical_jacobian(void)
x[i]->moles -= d2;
break;
}
if (use.Get_surface_ptr() != NULL)
{
base_surface = *use.Get_surface_ptr();
}
if (use.Get_gas_phase_ptr() != NULL)
{
*use.Get_gas_phase_ptr() = base_gas_phase;

View File

@ -1951,6 +1951,7 @@ jacobian_pz(void)
std::vector<class phase*> phase_ptrs;
std::vector<class phase> base_phases;
cxxGasPhase base_gas_phase;
cxxSurface base_surface;
LDBLE d, d1, d2;
int i, j;
Restart:
@ -1962,6 +1963,10 @@ Restart:
}
mb_sums();
residuals();
if (use.Get_surface_ptr() != NULL)
{
base_surface = *use.Get_surface_ptr();
}
if (use.Get_gas_phase_ptr() != NULL)
{
cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
@ -2141,6 +2146,10 @@ Restart:
reset();
break;
}
if (use.Get_surface_ptr() != NULL)
{
base_surface = *use.Get_surface_ptr();
}
if (use.Get_gas_phase_ptr() != NULL)
{
*use.Get_gas_phase_ptr() = base_gas_phase;
@ -2155,15 +2164,15 @@ Restart:
mb_sums();
residuals();
}
//for (i = 0; i < count_unknowns; i++)
//{
// //Debugging
// if (fabs(2.0 * (residual[i] - base[i]) / (residual[i] + base[i])) > 1e-2 &&
// fabs(residual[i]) + fabs(base[i]) > 1e-8)
// {
// std::cerr << i << ": " << x[i]->description << " " << residual[i] << " " << base[i] << std::endl;
// }
//}
for (i = 0; i < count_unknowns; i++)
{
//Debugging
if (fabs(2.0 * (residual[i] - base[i]) / (residual[i] + base[i])) > 1e-2 &&
fabs(residual[i]) + fabs(base[i]) > 1e-6)
{
std::cerr << i << ": " << x[i]->description << " " << residual[i] << " " << base[i] << std::endl;
}
}
base.clear();
calculating_deriv = 0;
return OK;

10
sit.cpp
View File

@ -967,7 +967,12 @@ jacobian_sit(void)
std::vector<class phase> base_phases;
double base_mass_water_bulk_x = 0, base_moles_h2o = 0;
cxxGasPhase base_gas_phase;
cxxSurface base_surface;
Restart:
if (use.Get_surface_ptr() != NULL)
{
base_surface = *use.Get_surface_ptr();
}
if (use.Get_gas_phase_ptr() != NULL)
{
cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
@ -1120,7 +1125,10 @@ Restart:
x[i]->moles -= d2;
break;
}
if (use.Get_surface_ptr() != NULL)
{
base_surface = *use.Get_surface_ptr();
}
if (use.Get_gas_phase_ptr() != NULL)
{
*use.Get_gas_phase_ptr() = base_gas_phase;