diff --git a/FEM/FEMEnums.h b/FEM/FEMEnums.h index 042823155..46b78b2a4 100644 --- a/FEM/FEMEnums.h +++ b/FEM/FEMEnums.h @@ -283,6 +283,14 @@ enum SolidReactiveSystem SolidReactiveSystem convertSolidReactiveSystem(const std::string& reactive_string); std::string convertSolidReactiveSystemToString(SolidReactiveSystem reactive_system); +enum InitDataReadWriteType +{ + NO_IO = 0, + READ, + WRITE, + READ_WRITE, +}; + } // end namespace FiniteElement struct TimType diff --git a/FEM/fem_ele.cpp b/FEM/fem_ele.cpp index 9bd0edab8..ab504b702 100644 --- a/FEM/fem_ele.cpp +++ b/FEM/fem_ele.cpp @@ -611,7 +611,8 @@ void CElement::computeJacobian(const int gp, const int order, const bool need_in DetJac = _Jacobian[0] * _Jacobian[3] - _Jacobian[1] * _Jacobian[2]; if (fabs(DetJac) < MKleinsteZahl) { - std::cout << "\n*** Jacobian: Det == 0 " << DetJac << "\n"; + std::cout << "\n*** Jacobian: Det == 0 " << DetJac + << " in element " << MeshElement->GetIndex() << "\n"; abort(); } invJacobian[0] = _Jacobian[3]; @@ -662,7 +663,8 @@ void CElement::computeJacobian(const int gp, const int order, const bool need_in if (fabs(DetJac) < MKleinsteZahl) { - std::cout << "\n*** Jacobian: DetJac == 0 " << DetJac << "\n"; + std::cout << "\n*** Jacobian: DetJac == 0 " << DetJac + << " in element " << MeshElement->GetIndex() << "\n"; abort(); } invJacobian[0] = _Jacobian[4] * _Jacobian[8] - _Jacobian[7] * _Jacobian[5]; diff --git a/FEM/fem_ele_std.cpp b/FEM/fem_ele_std.cpp index 09a166ec3..3fbb5ea49 100644 --- a/FEM/fem_ele_std.cpp +++ b/FEM/fem_ele_std.cpp @@ -1522,100 +1522,48 @@ void CFiniteElementStd::CalNodalEnthalpy() **************************************************************************/ double CFiniteElementStd::CalCoefMass() { - int Index = MeshElement->GetIndex(); + const int Index = MeshElement->GetIndex(); double val = 0.0; - double humi = 1.0; - double rhov = 0.0; - double arg[2]; - double biot_val, poro_val = 0.0, rho_val = 0.0, Se; - int tr_phase = 0; // SB, BG - double saturation = 0.0; // SB, BG - CompProperties* m_cp = NULL; - double drho_dp_rho = 0.0; // drho/dp * 1/rho - switch (PcsType) { default: std::cout << "Fatal error in CalCoefMass: No valid PCS type" << "\n"; + exit(EXIT_FAILURE); break; case EPT_LIQUID_FLOW: // Liquid flow - // Is this really needed? - val = MediaProp->StorageFunction(Index, unit, pcs->m_num->ls_theta); - - // get drho/dp/rho from material model or direct input - if (FluidProp->compressibility_model_pressure > 0) { - rho_val = FluidProp->Density(); - arg[0] = interpolate(NodalVal1); // p - arg[1] = interpolate(NodalValC1); // T - drho_dp_rho = FluidProp->drhodP(arg) / rho_val; - } - else - drho_dp_rho = FluidProp->drho_dp; - - // JT 2010, needed storage term and fluid compressibility... - // We derive here the storage at constant strain, or the inverse of Biot's "M" coefficient - // Assumptions are the most general possible:: Invarience under "pi" (Detournay & Cheng) loading. - // Se = 1/M = poro/Kf + (alpha-poro)/Ks :: Cf = 1/Kf = 1/rho * drho/dp :: alpha = 1 - K/Ks - // Second term (of Se) below vanishes for incompressible grains - // WW if(D_Flag > 0 && rho_val > MKleinsteZahl) - if (dm_pcs && MediaProp->storage_model == 7) // Add MediaProp->storage_model. 29.09.2011. WW - { - biot_val = SolidProp->biot_const; - poro_val = MediaProp->Porosity(Index, pcs->m_num->ls_theta); - val = 0.; // WX:04.2013 - - // WW if(SolidProp->K == 0) //WX: if HM Partitioned, K still 0 here - if (fabs(SolidProp->K) < DBL_MIN) // WW 29.09.2011 + val = MediaProp->StorageFunction(Index, unit, pcs->m_num->ls_theta); + double drho_dp_rho = 0.0; + // get drho/dp/rho from material model or direct input + if (FluidProp->compressibility_model_pressure > 0) { - if (SolidProp->Youngs_mode < 10 || SolidProp->Youngs_mode > 13) // JM,WX: 2013 - SolidProp->K = SolidProp->E / 3 / (1 - 2 * SolidProp->PoissonRatio); - else - { - double E_av; // average Youngs modulus - double nu_av; // average Poisson ratio - double nu_ai; // Poisson ratio perpendicular to the plane of isotropie, due to strain in the - // plane of isotropie - double nu_ia; // Poisson ratio in the plane of isotropie, due to strain perpendicular to the - // plane of isotropie - double nu_i; // Poisson ratio in the plane of isotropy - - E_av = 2. / 3. * (*SolidProp->data_Youngs)(0) + 1. / 3. * (*SolidProp->data_Youngs)(1); - - nu_ia = (*SolidProp->data_Youngs)(2); - nu_ai = nu_ia * (*SolidProp->data_Youngs)(1) - / (*SolidProp->data_Youngs)(0); // nu_ai=nu_ia*Ea/Ei - - nu_i = SolidProp->Poisson_Ratio(); - // 12 13 21 23 31 32 - // ai ai ia ii ia ii - nu_av = 1. / 3. * (nu_ai + nu_ia + nu_i); - - SolidProp->K = E_av / 3 / (1 - 2 * nu_av); - } + const double rho_val = FluidProp->Density(); + double arg[2]; + arg[0] = interpolate(NodalVal1); // p + arg[1] = interpolate(NodalValC1); // T + drho_dp_rho = FluidProp->drhodP(arg) / rho_val; + } + else + { + drho_dp_rho = FluidProp->drho_dp; } - val += poro_val * drho_dp_rho + (biot_val - poro_val) * (1.0 - biot_val) / SolidProp->K; - // Will handle the dual porosity version later... - } - else - { - poro_val = MediaProp->Porosity(Index, pcs->m_num->ls_theta); + const double poro_val = MediaProp->Porosity(Index, pcs->m_num->ls_theta); val += poro_val * drho_dp_rho; - } - // AS,WX: 08.2012 storage function eff stress - if (MediaProp->storage_effstress_model > 0) - { - double storage_effstress = 1.; - CFiniteElementStd* h_fem; - h_fem = this; - storage_effstress = MediaProp->StorageFunctionEffStress(Index, nnodes, h_fem); - val *= storage_effstress; - } + // AS,WX: 08.2012 storage function eff stress + if (MediaProp->storage_effstress_model > 0) + { + double storage_effstress = 1.; + CFiniteElementStd* h_fem; + h_fem = this; + storage_effstress = MediaProp->StorageFunctionEffStress(Index, nnodes, h_fem); + val *= storage_effstress; + } - //val /= time_unit_factor; + //val /= time_unit_factor; + } break; case EPT_UNCONFINED_FLOW: // Unconfined flow break; @@ -1640,12 +1588,12 @@ double CFiniteElementStd::CalCoefMass() // JT 2010, generalized poroelastic storage. See single phase version in case "L". // Se = 1/M = poro/Kf + (alpha-poro)/Ks - rho_val = FluidProp->Density(); + const double rho_val = FluidProp->Density(); if (D_Flag > 0 && rho_val > MKleinsteZahl) { - biot_val = SolidProp->biot_const; - poro_val = MediaProp->Porosity(Index, pcs->m_num->ls_theta); - Se = poro_val * FluidProp->drho_dp + (biot_val - poro_val) * (1.0 - biot_val) / SolidProp->K; + const double biot_val = SolidProp->biot_const; + const double poro_val = MediaProp->Porosity(Index, pcs->m_num->ls_theta); + const double Se = poro_val * FluidProp->drho_dp + (biot_val - poro_val) * (1.0 - biot_val) / SolidProp->K; // The poroelastic portion val += Se * MMax(0., Sw); } @@ -1660,6 +1608,7 @@ double CFiniteElementStd::CalCoefMass() // dSedPc always return positive numbers for default case // However, the value should be negative analytically. double dSwdPc = m_mmp->PressureSaturationDependency(Sw, true); + const double poro_val = MediaProp->Porosity(Index, pcs->m_num->ls_theta); val += poro_val * dSwdPc; } } @@ -1682,65 +1631,71 @@ double CFiniteElementStd::CalCoefMass() break; //.................................................................... case EPT_MASS_TRANSPORT: // Mass transport //SB4200 - // Porosity - val = MediaProp->Porosity(Index, pcs->m_num->ls_theta); - // SB Transport in both phases - tr_phase = cp_vec[this->pcs->pcs_component_number]->transport_phase; - // Multi phase transport of components - saturation = PCSGetEleMeanNodeSecondary_2(Index, pcs->flow_pcs_type, "SATURATION1", 1); - if (tr_phase == 0) // Water phase - val *= saturation; - else if (tr_phase == 10) // non wetting phase - val *= (1.0 - saturation); - m_cp = cp_vec[pcs->pcs_component_number]; - // Retardation Factor - val *= m_cp->CalcElementRetardationFactorNew(Index, unit, pcs); + { + // Porosity + val = MediaProp->Porosity(Index, pcs->m_num->ls_theta); + // SB Transport in both phases + const double tr_phase = cp_vec[this->pcs->pcs_component_number]->transport_phase; + // Multi phase transport of components + const double saturation = PCSGetEleMeanNodeSecondary_2(Index, pcs->flow_pcs_type, "SATURATION1", 1); + if (tr_phase == 0) // Water phase + val *= saturation; + else if (tr_phase == 10) // non wetting phase + val *= (1.0 - saturation); + CompProperties* const m_cp = cp_vec[pcs->pcs_component_number]; + // Retardation Factor + val *= m_cp->CalcElementRetardationFactorNew(Index, unit, pcs); + } break; case EPT_OVERLAND_FLOW: // Liquid flow val = 1.0; break; case EPT_RICHARDS_FLOW: // Richards - Sw = 1.0; - dSdp = 0.; - PG = interpolate(NodalVal1); // 12.02.2007. Important! WW - - if (PG < 0.0) // JM skip to call these two functions in saturated case { - Sw = MediaProp->SaturationCapillaryPressureFunction(-PG); - dSdp = -MediaProp->PressureSaturationDependency(Sw, true); // JT: dSdp now returns actual sign (i.e. <0) - } + Sw = 1.0; + dSdp = 0.; + PG = interpolate(NodalVal1); // 12.02.2007. Important! WW - poro = MediaProp->Porosity(Index, pcs->m_num->ls_theta); - rhow = FluidProp->Density(); + if (PG < 0.0) // JM skip to call these two functions in saturated case + { + Sw = MediaProp->SaturationCapillaryPressureFunction(-PG); + dSdp = -MediaProp->PressureSaturationDependency(Sw, true); // JT: dSdp now returns actual sign (i.e. <0) + } - if (FluidProp->compressibility_model_pressure > 0) - { // drho_dp from rho-p-T relation - arg[0] = PG; - arg[1] = interpolate(NodalValC1); // T - drho_dp_rho = FluidProp->drhodP(arg) / rhow; - } - else - drho_dp_rho = FluidProp->drho_dp; + poro = MediaProp->Porosity(Index, pcs->m_num->ls_theta); + rhow = FluidProp->Density(); - // Storativity - val = MediaProp->StorageFunction(Index, unit, pcs->m_num->ls_theta) * Sw; + double drho_dp_rho = 0.0; + if (FluidProp->compressibility_model_pressure > 0) + { // drho_dp from rho-p-T relation + double arg[2]; + arg[0] = PG; + arg[1] = interpolate(NodalValC1); // T + drho_dp_rho = FluidProp->drhodP(arg) / rhow; + } + else + drho_dp_rho = FluidProp->drho_dp; - // Fluid compressibility - if (rhow > 0.0) - val += poro * Sw * drho_dp_rho; - // Capillarity - if (PG < 0.0) // dSdp gives always a value>0, even if p>0! - val += poro * dSdp; - // WW - if (MediaProp->heat_diffusion_model == 1) - { - // PG = fabs(interpolate(NodalVal1)); - TG = interpolate(NodalValC) + PhysicalConstant::CelsiusZeroInKelvin; - humi = exp(PG / (SpecificGasConstant::WaterVapour * TG * rhow)); - rhov = humi * FluidProp->vaporDensity(TG); - // - val -= poro * rhov * dSdp / rhow; - val += (1.0 - Sw) * poro * rhov / (rhow * rhow * SpecificGasConstant::WaterVapour * TG); + // Storativity + val = MediaProp->StorageFunction(Index, unit, pcs->m_num->ls_theta) * Sw; + + // Fluid compressibility + if (rhow > 0.0) + val += poro * Sw * drho_dp_rho; + // Capillarity + if (PG < 0.0) // dSdp gives always a value>0, even if p>0! + val += poro * dSdp; + // WW + if (MediaProp->heat_diffusion_model == 1) + { + // PG = fabs(interpolate(NodalVal1)); + TG = interpolate(NodalValC) + PhysicalConstant::CelsiusZeroInKelvin; + const double humi = exp(PG / (SpecificGasConstant::WaterVapour * TG * rhow)); + const double rhov = humi * FluidProp->vaporDensity(TG); + // + val -= poro * rhov * dSdp / rhow; + val += (1.0 - Sw) * poro * rhov / (rhow * rhow * SpecificGasConstant::WaterVapour * TG); + } } break; case EPT_FLUID_MOMENTUM: // Fluid Momentum diff --git a/FEM/pcs_dm.cpp b/FEM/pcs_dm.cpp index 4a09638a0..726659e25 100644 --- a/FEM/pcs_dm.cpp +++ b/FEM/pcs_dm.cpp @@ -9,8 +9,6 @@ #include "pcs_dm.h" -#include "makros.h" - #include #include #include @@ -18,6 +16,9 @@ #include #include +#include "display.h" +#include "makros.h" + #include "StringTools.h" #include "FEMEnums.h" @@ -89,7 +90,7 @@ using Math_Group::Matrix; namespace process { CRFProcessDeformation::CRFProcessDeformation() - : CRFProcess(), fem_dm(NULL), ARRAY(NULL), counter(0), InitialNorm(0.0), idata_type(none), + : CRFProcess(), fem_dm(NULL), ARRAY(NULL), counter(0), InitialNorm(0.0), _has_initial_stress_data(false), error_k0(1.0e10) { @@ -102,7 +103,8 @@ CRFProcessDeformation::~CRFProcessDeformation() if (myrank == 0) #endif WriteGaussPointStress(last_step); - if (type == 41 && (idata_type == write_all_binary || idata_type == read_write)) + if (type == 41 && (_init_domain_data_type == FiniteElement::WRITE + || _init_domain_data_type == FiniteElement::READ_WRITE)) { // mono-deformation-liquid WriteSolution(); @@ -152,16 +154,6 @@ CRFProcessDeformation::~CRFProcessDeformation() **************************************************************************/ void CRFProcessDeformation::Initialization() { - //-- NW 25.10.2011 - // this section has to be executed at latest before calling InitGauss() - // Control for reading and writing solution - if (reload == 1) - idata_type = write_all_binary; - if (reload == 2) - idata_type = read_all_binary; - if (reload == 3) - idata_type = read_write; - // Local assembliers // An instaniate of CFiniteElementVec int i, Axisymm = 1; // ani-axisymmetry @@ -1214,7 +1206,8 @@ void CRFProcessDeformation::InitGauss(void) } } // Reload the stress results of the previous simulation - if (idata_type == read_all_binary || idata_type == read_write) + if ( _init_domain_data_type == FiniteElement::READ + || _init_domain_data_type == FiniteElement::READ_WRITE) { ReadGaussPointStress(); if (type == 41) // mono-deformation-liquid @@ -2986,7 +2979,8 @@ void CRFProcessDeformation::UpdateStress() **************************************************************************/ void CRFProcessDeformation::WriteGaussPointStress(const bool last_step) { - if (!(idata_type == write_all_binary || idata_type == read_write)) + if (!(_init_domain_data_type == FiniteElement::WRITE + || _init_domain_data_type == FiniteElement::READ_WRITE)) return; if ((aktueller_zeitschritt % nwrite_restart) > 0 && (!last_step)) @@ -2997,7 +2991,7 @@ void CRFProcessDeformation::WriteGaussPointStress(const bool last_step) #else const std::string StressFileName = FileName + ".sts"; #endif - + ScreenMessage("-> Write initial stress \n"); fstream file_stress(StressFileName.data(), ios::binary | ios::out | ios::trunc); ElementValue_DM* eleV_DM = NULL; @@ -3049,6 +3043,7 @@ void CRFProcessDeformation::ReadGaussPointStress() const std::string StressFileName = FileName + ".sts"; #endif + ScreenMessage("-> Read initial stress \n"); fstream file_stress(StressFileName.data(), ios::binary | ios::in); ElementValue_DM* eleV_DM = NULL; // diff --git a/FEM/pcs_dm.h b/FEM/pcs_dm.h index e8f13a1e6..1d6dedb36 100644 --- a/FEM/pcs_dm.h +++ b/FEM/pcs_dm.h @@ -18,6 +18,7 @@ #include +#include "FEMEnums.h" #include "rf_pcs.h" // Strong discontinuity @@ -43,22 +44,6 @@ class CPARDomain; namespace process { -enum InitDataReadWriteType -{ - none, - read_write, - read_all_binary, - write_all_binary, - read_all_asci, - write_all_asci, - read_stress_binary, - write_stress_binary, - read_displacement, - write_displacement, - read_pressure, - write_pressure -}; - // Elasto-plastic Deformation class CRFProcessDeformation : public CRFProcess { @@ -138,14 +123,12 @@ class CRFProcessDeformation : public CRFProcess double InitialNormU; double InitialNormU0; - InitDataReadWriteType idata_type; - bool _has_initial_stress_data; // double error_k0; #if !defined(USE_PETSC) // && !defined(other parallel libs)//03.3012. WW - // Domain decompisition + // Domain decomposition void DomainAssembly(CPARDomain* m_dom); #endif diff --git a/FEM/problem.cpp b/FEM/problem.cpp index 1a7e03e09..0dcdb195d 100644 --- a/FEM/problem.cpp +++ b/FEM/problem.cpp @@ -1741,7 +1741,8 @@ void Problem::PostCouplingLoop() // the first pcs process is the flow process...if reload not defined for every process, // restarting with gems will not work in any case - if ((m_pcs->reload == 1 || m_pcs->reload == 3) + if (( m_pcs->GetRestartFlag() == FiniteElement::WRITE + || m_pcs->GetRestartFlag() == FiniteElement::READ_WRITE) && !((aktueller_zeitschritt % m_pcs->nwrite_restart) > 0)) m_vec_GEM->WriteReloadGem(); diff --git a/FEM/rf_REACT_GEM.cpp b/FEM/rf_REACT_GEM.cpp index 7ab8192a0..27077ce1d 100644 --- a/FEM/rf_REACT_GEM.cpp +++ b/FEM/rf_REACT_GEM.cpp @@ -794,7 +794,9 @@ short REACT_GEM::Init_RUN(string Project_path) // TODO for petsc we may need to synchronize saturations ...have to check! } - if ((m_flow_pcs->GetRestartFlag() >= 2)) // get the restart data specific for gems + // get the restart data specific for gems + if ( m_flow_pcs->GetRestartFlag() == FiniteElement::READ + || m_flow_pcs->GetRestartFlag() == FiniteElement::READ_WRITE) { if (!ReadReloadGem()) { @@ -964,7 +966,9 @@ short REACT_GEM::GetInitialReactInfoFromMassTransport(int timelevel) // get pressure from MT m_P[node_i] = REACT_GEM::GetPressureValue_MT(node_i, timelevel); // get Independent and dependent Component value from MT - if ((flag_transport_b == 1) && (m_flow_pcs->GetRestartFlag() < 2)) + if ((flag_transport_b == 1) && + ( m_flow_pcs->GetRestartFlag() == FiniteElement::WRITE + || m_flow_pcs->GetRestartFlag() == FiniteElement::NO_IO)) REACT_GEM::GetBValue_MT(node_i, timelevel, m_bIC + node_i * nIC); // do this not for restart...overwrites values!!! } @@ -3197,9 +3201,9 @@ int REACT_GEM::CalcLimitsInitial(long in, TNode* m_Node) m_dll[in * nDC + j] = 0.0; // set to zero m_dul[in * nDC + j] = 1.0e+10; // very high number } - - if ((m_flow_pcs->GetRestartFlag() - >= 2)) // we test if restart flag is set....if not this will not work, as x_dc might be not correct + // we test if restart flag is set....if not this will not work, as x_dc might be not correct + if ( m_flow_pcs->GetRestartFlag() == FiniteElement::READ + || m_flow_pcs->GetRestartFlag() == FiniteElement::READ_WRITE) { for (ii = 0; ii < (int)m_kin.size(); ii++) { @@ -4201,16 +4205,16 @@ void REACT_GEM::gems_worker(int tid, string m_Project_path) // rwmutex.lock(); // cout << "GEMS3K MPI Processe / Thread: " << myrank << " " << tid << " in " << in << "\n"; // rwmutex.unlock(); - - if ((m_flow_pcs->GetRestartFlag() - >= 2)) // everything is stored in concentrations for restart ...moved it to here from init_gems - + // everything is stored in concentrations for restart ...moved it to here from init_gems + if ( m_flow_pcs->GetRestartFlag() == FiniteElement::READ + || m_flow_pcs->GetRestartFlag() == FiniteElement::READ_WRITE) // Convert from concentration REACT_GEM::ConcentrationToMass(in, 1); // I believe this is save for MPI // this we have already // now we calculate kinetic constraints for GEMS! - if (!(m_flow_pcs->GetRestartFlag() >= 2)) + if ( m_flow_pcs->GetRestartFlag() == FiniteElement::WRITE + || m_flow_pcs->GetRestartFlag() == FiniteElement::NO_IO) REACT_GEM::CalcLimitsInitial( in, t_Node); // kg44 16.05.2013 new version, restart files contain upper and lower limits // Manipulate some kinetic contstraints for special initial conditions @@ -4292,8 +4296,9 @@ void REACT_GEM::gems_worker(int tid, string m_Project_path) // calculate density of fluid phase, which is normally the first phase m_fluid_density[in] = m_mPS[in * nPS + 0] / m_vPS[in * nPS + 0]; - if ((m_flow_pcs->GetRestartFlag() < 2)) // we do not need the second pass for complete restart - + // we do not need the second pass for complete restart + if ( m_flow_pcs->GetRestartFlag() == FiniteElement::WRITE + || m_flow_pcs->GetRestartFlag() == FiniteElement::NO_IO) // scale data so that second pass gives the normalized volume of 1m^3 if (m_Vs[in] >= 0.0) { // this should be not done for restart, decoupled porosity runs do not change volumes and the other runs @@ -4372,7 +4377,8 @@ void REACT_GEM::gems_worker(int tid, string m_Project_path) // } // end if for restart // calculate the chemical porosity - if (m_flow_pcs->GetRestartFlag() < 2) + if ( m_flow_pcs->GetRestartFlag() == FiniteElement::WRITE + || m_flow_pcs->GetRestartFlag() == FiniteElement::NO_IO) REACT_GEM::CalcPorosity(in, t_Node); // during init it should be always done, except for restart !!! REACT_GEM::CalcReactionRate( diff --git a/FEM/rf_mmp_new.cpp b/FEM/rf_mmp_new.cpp index 2b59452b8..e933f6450 100644 --- a/FEM/rf_mmp_new.cpp +++ b/FEM/rf_mmp_new.cpp @@ -774,12 +774,13 @@ std::ios::pos_type CMediumProperties::Read(std::ifstream* mmp_file) << "\n"; break; case 7: // RW/WW - in >> storage_model_values[0]; // Biot's alpha - in >> storage_model_values[1]; // Skempton's B coefficient - in >> storage_model_values[2]; // macroscopic drained bulk modulus - double val_l = storage_model_values[0] * (1. - storage_model_values[0] * storage_model_values[1]) - / storage_model_values[1] / storage_model_values[2]; - storage_model_values[1] = val_l; + { + if(!PCSGet("DEFORMATION")) + { + ScreenMessage("Error: Porosity model 7 must be combined with deformation process"); + exit(EXIT_FAILURE); + } + } break; } in.clear(); @@ -7351,8 +7352,19 @@ double CMediumProperties::StorageFunction(long index, double* gp, double theta) */ break; case 7: // poroelasticity RW - storage = storage_model_values[1]; - break; + { + // Moved the following comment from double CFiniteElementStd::CalCoefMass() to here by WW + // JT 2010, needed storage term and fluid compressibility... + // We derive here the storage at constant strain, or the inverse of Biot's "M" coefficient + // Assumptions are the most general possible:: Invarience under "pi" (Detournay & Cheng) loading. + // Se = 1/M = poro/Kf + (alpha-poro)/Ks :: Cf = 1/Kf = 1/rho * drho/dp :: alpha = 1 - K/Ks + // Second term (of Se) below vanishes for incompressible grains + + SolidProp::CSolidProperties* solid_prop = msp_vector[Fem_Ele_Std->GetMeshElement()->GetPatchIndex()]; + const double biots_constant = solid_prop->getBiotsConstant(); + const double porosity = Porosity(index, theta); + return (biots_constant - porosity) * (1.0 - biots_constant) / solid_prop->getBulkModulus(); + } case 10: if (permeability_saturation_model[0] == 10) // MW storage = porosity_model_values[0] / (gravity_constant * gravity_constant * mfp_vector[0]->Density()); diff --git a/FEM/rf_msp_new.cpp b/FEM/rf_msp_new.cpp index e3581a339..819e8a161 100644 --- a/FEM/rf_msp_new.cpp +++ b/FEM/rf_msp_new.cpp @@ -8180,6 +8180,8 @@ double CSolidProperties::E_Function(int dim, const ElementValue_DM* ele_val, int return_value = GetCurveValue((int)E_Function_Model_Value[0], 0, prin_str[0], &valid); break; } + case 3: + return GetCurveValue((int)E_Function_Model_Value[0], 0, aktuelle_zeit, &valid); default: return_value = 1.; break; @@ -8213,6 +8215,38 @@ double CSolidProperties::getBishopCoefficient(const double effectiveS, const dou } return p; } + +double CSolidProperties::getBulkModulus() const +{ + if (K > DBL_MIN) + return K; + + if (Youngs_mode < 10 || Youngs_mode > 13) + return E / 3 / (1 - 2 * PoissonRatio); + + // average Youngs modulus + double const E_av = 2. / 3. * (*data_Youngs)(0) + 1. / 3. * (*data_Youngs)(1); + + // Poisson ratio perpendicular to the plane of isotropie, due to strain in the + // plane of isotropie + const double nu_ia = (*data_Youngs)(2); + + // Poisson ratio in the plane of isotropie, due to strain perpendicular to the + // plane of isotropie + const double nu_ai = nu_ia * (*data_Youngs)(1) + / (*data_Youngs)(0); // nu_ai=nu_ia*Ea/Ei + + // Poisson ratio in the plane of isotropy + const double nu_i = Poisson_Ratio(); + + // average Poisson ratio + // 12 13 21 23 31 32 + // ai ai ia ii ia ii + const double nu_av = 1. / 3. * (nu_ai + nu_ia + nu_i); + + return E_av / 3 / (1 - 2 * nu_av); +} + } // end namespace ///////////////////////////////////////////////////////////////////////////// diff --git a/FEM/rf_msp_new.h b/FEM/rf_msp_new.h index ccf449203..fd3d92b6f 100644 --- a/FEM/rf_msp_new.h +++ b/FEM/rf_msp_new.h @@ -163,6 +163,10 @@ class CSolidProperties std::vector& eps_K_curr, std::vector& eps_M_curr, std::vector& eps_pl_curr, double& e_pl_v, double& e_pl_eff, double& lam, Math_Group::Matrix& Consistent_Tangent, double Temperature, double& local_res); + + double getBulkModulus() const; + double getBiotsConstant() const { return biot_const; } + private: // CMCD FiniteElement::CFiniteElementStd* Fem_Ele_Std; diff --git a/FEM/rf_pcs.cpp b/FEM/rf_pcs.cpp index e68d56339..51d7804ee 100644 --- a/FEM/rf_pcs.cpp +++ b/FEM/rf_pcs.cpp @@ -235,7 +235,7 @@ T* resize(T* array, size_t old_size, size_t new_size); CRFProcess::CRFProcess(void) : _problem(NULL), p_var_index(NULL), num_nodes_p_var(NULL), fem(NULL), Memory_Type(0), Write_Matrix(false), matrix_file(NULL), - WriteSourceNBC_RHS(0), + WriteSourceNBC_RHS(0), _init_domain_data_type(FiniteElement::NO_IO), #ifdef JFNK_H2M JFNK_precond(false), norm_u_JFNK(NULL), array_u_JFNK(NULL), array_Fu_JFNK(NULL), #endif @@ -322,8 +322,6 @@ CRFProcess::CRFProcess(void) selected = true; // OK // MSH OK m_msh = NULL; - // Reload solutions - reload = -1; nwrite_restart = 1; // kg44 write every timestep is default pcs_nval_data = NULL; pcs_eval_data = NULL; @@ -544,6 +542,13 @@ CRFProcess::~CRFProcess(void) delete eqs_new; eqs_new = NULL; #endif + if (_init_domain_data_type == FiniteElement::READ) + { + ScreenMessage("!!!Waning: Reading and writing solutions (reload=3)\ + are selected in this computation.\n"); + ScreenMessage("\t Please make sure that the original initial\ + data files are used.\n"); + } } /************************************************************************** @@ -870,7 +875,9 @@ void CRFProcess::Create() } } // - if (reload >= 2 && type != 4 && type / 10 != 4) // Modified at 03.08.2010. WW + if (( _init_domain_data_type == FiniteElement::READ + || _init_domain_data_type == FiniteElement::READ_WRITE) + && type != 4 && type / 10 != 4) // Modified at 03.08.2010. WW { // PCH cout << "Reloading the primary variables... " @@ -878,7 +885,8 @@ void CRFProcess::Create() ReadSolution(); // WW } - if (reload < 2) // PCH: If reload is set, no need to have ICs + if ( _init_domain_data_type == FiniteElement::NO_IO + || _init_domain_data_type == FiniteElement::WRITE ) // PCH: If reload is set, no need to have ICs { // IC cout << "->Assign IC" << '\n'; @@ -886,8 +894,8 @@ void CRFProcess::Create() } else // Bypassing IC - cout << "RELOAD is set to be " << reload << ". So, bypassing IC's" - << "\n"; + cout << "RELOAD is set to be " << _init_domain_data_type + << ". So, bypassing IC's" << "\n"; if (pcs_type_name_vector.size() && pcs_type_name_vector[0].find("DYNAMIC") != string::npos) // WW setIC_danymic_problems(); @@ -1235,7 +1243,8 @@ void CRFProcess::Write_Processed_BC() **************************************************************************/ void CRFProcess::WriteSolution() { - if (reload == 2 || reload <= 0) + if ( _init_domain_data_type == FiniteElement::NO_IO + || _init_domain_data_type == FiniteElement::READ) return; // kg44 write out only between nwrite_restart timesteps if ((aktueller_zeitschritt % nwrite_restart) > 0) @@ -1719,7 +1728,7 @@ CRFProcess* CRFProcess::CopyPCStoDM_PCS() dm_pcs->num_type_name = num_type_name; dm_pcs->Memory_Type = Memory_Type; dm_pcs->NumDeactivated_SubDomains = NumDeactivated_SubDomains; - dm_pcs->reload = reload; + dm_pcs->_init_domain_data_type = _init_domain_data_type; dm_pcs->nwrite_restart = nwrite_restart; dm_pcs->isPCSDeformation = true; dm_pcs->isPCSFlow = this->isPCSFlow; // JT @@ -2003,7 +2012,14 @@ std::ios::pos_type CRFProcess::Read(std::ifstream* pcs_file) // subkeyword found if (line_string.find("$RELOAD") != string::npos) { + int reload; *pcs_file >> reload; // WW + if (reload == 1) + _init_domain_data_type = FiniteElement::WRITE; + if (reload == 2) + _init_domain_data_type = FiniteElement::READ; + if (reload == 3) + _init_domain_data_type = FiniteElement::READ_WRITE; if (reload == 1 || reload == 3) *pcs_file >> nwrite_restart; // kg44 read number of timesteps between writing restart files pcs_file->ignore(MAX_ZEILE, '\n'); @@ -12530,7 +12546,9 @@ bool CRFProcess::NODRelations() // IC cout << "->Assign IC" << '\n'; - if (reload == 2 && type != 4 && type != 41) + if ( ( _init_domain_data_type == FiniteElement::READ + || _init_domain_data_type == FiniteElement::READ_WRITE) + && type != 4 && type != 41) ReadSolution(); // WW SetIC(); diff --git a/FEM/rf_pcs.h b/FEM/rf_pcs.h index a071895f1..d2eb7c4f4 100644 --- a/FEM/rf_pcs.h +++ b/FEM/rf_pcs.h @@ -320,11 +320,8 @@ class CRFProcess : public ProcessInfo int WriteSourceNBC_RHS; int WriteProcessed_BC; // Write the current solutions/Read the previous solutions WW - // -1, default. Do nothing - // 1. Write - // 2. Read - // 3 read and write - int reload; + FiniteElement::InitDataReadWriteType _init_domain_data_type; + long nwrite_restart; void WriteRHS_of_ST_NeumannBC(); void ReadRHS_of_ST_NeumannBC(); @@ -746,7 +743,8 @@ class CRFProcess : public ProcessInfo // WW void CheckSTGroup(); //OK #ifdef GEM_REACT void IncorporateSourceTerms_GEMS(void); // HS: dC/dt from GEMS chemical solver. - int GetRestartFlag() const { return reload; } + FiniteElement::InitDataReadWriteType GetRestartFlag() const + { return _init_domain_data_type; } #endif // BC void IncorporateBoundaryConditions(const int rank = -1);