Hi Mohammad,
I am surprised that you see different behaviour for porosity and permeability. At first glance I would have at least expected the same behaviour for the two.
Besides Timos comment: your code snippet from the main file shows that you register the permeability field, then tell your problem to update the field, and write out the initial conditions. During your simulation, are you updating the field after each time step (i.e. calling problem->updateVTKOutput())? I would guess that that's necessary, but I didn't see that in your explanations...
Cheers, Dennis On 12.10.22 12:40, Mohammad Hodroj (Student) wrote:
Dear DuMux Community, I hope this email finds you well. I am writing to kindly ask for your support in this matter: *_Given:_* 1. I am using the porous medium flow model (injection problem) 2. The Problem is : 3. template<class FVGridGeometry, class Scalar> class InjectionSpatialParams : public FVSpatialParams<*FVGridGeometry*, Scalar, InjectionSpatialParams<FVGridGeometry, Scalar>> 4. 5. PermeabilityType = Scalar ** 6. I am using PorosityAtPos (const GlobalPos& globalPos) const ** 5. { if (isInCement_(globalPos)) { // certain conditions cementPorosity_ = equation (xxx) // at each globalpos[2] there is a new porosity value return cementPorosity; } else if (isInAquifer_(globalPos)) return aquiferPorosity_; else return overburdenPorosity_; } 5. I am using PermeabilityAtPos (const GlobalPos& globalPos) const 6. 7. 6. { Scalar cementPor_; Scalar cementK_; const Scalar dc_Porosity = 0.3; Scalar InitialCementK_ = 8.645e-13; if (isInCement_(globalPos)) {cementPor_ = porosityAtPos( globalPos); // calling porosity function cementK_ = *InitialCementK_* *(std::pow(((1-dc_Porosity)/(1- cementPor_)),2) *(std::pow((cementPor_/dc_Porosity), 3))); //std::cout<< " Porosity:"<< cementPor_ <<" Permeability:"<< cementK_ << std::endl;return cementK_; } else if (isInAquifer_(globalPos)) return aquiferK_; else return overburdenK_;} 6. I am using a function to return permeability for additional vtk output : 7. 8. const std::vector<Scalar>& getpermeability() 9. { return permeability; } 7. I am using a function that updates the permeability for additional vtk output void updateVtkOutput(const SolutionVector& curSol) { auto fvGeometry = localView(this->gridGeometry());for (const auto& element : elements(this->gridGeometry().gridView())){const auto elemSol = elementSolution(element, curSol, this->gridGeometry());fvGeometry.bindElement(element); for (auto&& scv : scvs(fvGeometry)) { VolumeVariables volVars; volVars.update(elemSol, *this, element, scv); const auto dofIdxGlobal = scv.dofIndex(); permeability_[dofIdxGlobal] = volVars.permeability(); } } } 8. I added the following to return the permeability for addtional vtk output in main.cc: vtkWriter.addField(problem->getpermeability(), "Permeability"); problem->updateVtkOutput(x); vtkWriter.write(0.0); 9. I am using also the following:unsigned int codim = GetPropType<TypeTag, Properties::GridGeometry>::discMethod == DiscretizationMethod::box? dim : 0;permeability_.resize(fvGridGeometry->gridView().size(codim)); *Issue* 1. I can see the different permeability values of the cement layer (which are dependent on porosity)while printing them in the terminal but not in vtu files and ParaView. 2. I can see in the vtu files only the initial value of permeability (constant) of the cement layer.* Note: **I can see the new porosity values of the cement layer in ParaView as well as upon printing them in the terminal (porosity changing with time)* So, how can I get the new values of permeability as appearing in the terminal in the vtu files ( similar to porosity)?Thanks for your help. Best Regards, Mohammad Hodroj _______________________________________________ DuMux mailing list [email protected] https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
_______________________________________________ DuMux mailing list [email protected] https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
