PorousFlow: masking areas for conductive flow/transport

[garciapintado]

Dear all,
I am using the PorousFlow module, and wondering if you could advise on a way to simulate that some blocks are completely impervious. That is, I would like thermal diffusion to act in all the domain, but to be able to constrain the blocks for Darcy flow, so that conductive heat and advective transport in only simulated in the permeable domain. To use a negligible porosity [let's say 1.0E-2] in the assumed impervious areas prevents convergence in my attempts.
Many thanks,
Javier GP

[GiudGiud]

Hello

You should be able to block-restrict some physics. I'm not sure how it's done for porous flow. What does your input file look like? Do you have a large [Kernels] section with each term in the equations or is there a [PorousFlow] block (action) that handles the equations ?

Guillaume

[garciapintado]

Hello Guillaume,
Many thanks for the reply [sorry, I was aiming to ask for a general Q&A Modules, not specifically for Navier-Stokes].
I am new to MOOSE, and trying to find my way through it.
The PorousFlow modules has several actions, which handle the equations.
I’m using the PorousFlowFullySaturated action. Here, the heat equation in this includes (apart from the mechanical advection) both heat conduction and heat advection by Darcy flow (each with its specific kernel).
Following your suggestion, if you could point to some example in MOOSE (not necessarily related to the PorousFlow module) where some of the physics are restricted to specific blocks this would be great help. Basically, I'd need some blocks to be only subject to head conduction.
Javier

[GiudGiud]

No worries.

About block restriction of physics, here's an example:
https://github.com/idaholab/moose/blob/06bc1f73f0d4dfb3d01ce9ba77d6bb97346e07eb/test/tests/kernels/block_kernel/block_kernel_test.i
the body force is restricted to only a single subdomain.
But it's not always that simple, sometimes you'll need to add boundary conditions to restrict physics to some subdomains. It's the case for finite volume navier stokes (not the same as porous flow as you know)

But please see @WilkAndy 's answer

[WilkAndy]

The easiest thing is to set permeability zero or very low in those impermeable blocks.

Almost definitely this will work well, but before doing this, just make sure your physics is going to work properly. For instance, it sounds like you're doing a thermal-hydraulic simulation. If you set permeability = 0 in some blocks, and then reduce/increase the temperature when using a high-precision equation of state for your fluid, then you might start to probe the limits of the EOS's validity.

[garciapintado]

Hello Andy,
Just missed your last message. Well, you see I’ve been doing simple TH synthetic tests.
I have not modified any C++ code in MOOSe by now [actually I only know the basic of C++]. I had planned that if these tests work fine (so we go ahead with MOOSE), I’d add a new material with the parametrization I’ve put in Matlab, where for a start I used your own PorousFlowPorosity formulation [the fluid, thermal and mechanical parts], and modified some minor things that did not seem to fit our problem. Basically, I just used the volumetric
deformation [the dynamical pressure term] for the porosity mechanical factor, and then I added a factor dependent of the strain rate [i.e. its 2nd invariant] for the permeability, after apply a Kozeny-Carman formulation [also borrowed from your own code]. This seems to do what we expected in general terms in our Matlab code. I could add the details, but it s a quite simple thing.

So, in our case, the mechanics would change porosity [and permeability], and this is key for our simulations to work and properly develop the fluid convection cells.

About the  “T>800degC + porosity idea”, you’ve got me. I need to look at that…

I am enthusiastic about the prospect of MOOSE, if each 10 kyr of hydrothermal simulations [i.e. a single timestep of the mechanical model] could go down to about 3 minutes computing time within a single computing node this would allow us to do the runs [including ensembles] we need! [each of these 10 kyr now takes us 30’ in Matlab with similar number of processors, which make the computing impractical].

I’m going to check what I have for the high-temperature porosity…

Javier
[needless to say that if you find the application interesting and we manage to put this into work in MOOSE, we would be more than happy if you did not mind to coauthor derived results]. I also find myself the application fascinating...

[garciapintado]

I have checked. In the simulation for the plot I sent for our real case [the matlab simulation]. I just had the trace of the volumetric strain contributing to the mechanical effect on the porosity:
volumetric_strain = -1 / (bulk_modulus * (Pd - Pd_ref)),
where Pd is the mechanical model pressure. The fluid and thermal effect as in PorousflowPorosity, and then the simplistic approach of setting a fixed positive lower threshold of the resulting porosity for vey low values [I’m aware this is not good convergence-wise in MOOSE], which I guess is what you were asking for. So, no great idea here. But as I've seen I believe we rarely hit these low values in our zone of interest.

Then, for the permeability, the Kozeny-Carman permeability, [after your kozeny_carman_phi0] and after this permeability (K) has been calculated, it is multiplied by a factor which result from a function [a Gaussian-like increase, as the classical Gaussian variograms] dependent on the ratio of the [IInd invariant] of the plastic strain rate to a nominal value [which is around the maximum value we normally have around big faults with highly localised deformation for this class of model].

All in all, the permeability becomes clearly higher on active faults in the brittle domain, and then decrease when the faults become inactive. Part of this idea is to emulate the fault remineralisation effect that we are lacking in our model. That is, we cannot use total accumulated strain to control permeability as we would be neglecting the remineralisation in abandoned faults which would seal them for our geological timescales.

[WilkAndy]

Reading the above, i believe that MOOSE is behaving OK with SimpleFluidProperties. Tell me if that's not true.

I'm sure that the zero derivatives you have in your extrapolation will negatively impact convergence. Even if the final converged solution is within reasonable bounds, the PETSc/MOOSE solver will always sample some crazy values of (P, T) during the convergence process. If you have zero derivative you're not giving PETSc/MOOSE any clues about how density and viscosity change with (P, T), so it thinks it can choose random values and that density and viscosity will remain constant, so it chooses those random values (like T = -1500.82) and things go crazy. (Actually, it may not be quite that silly, but that's the principle anyway.)

Do you think your trials indicate we simply must be more clever in the Tabulated IAPWS?

(Aside: i try to have about 10k DOFs per compute node.)

How exactly are you planning to feed the mechanical stuff into your MOOSE simulations? Perhaps as some files defining the spatial permeability and porosity?

After getting those trial simulations mentioned above working (i think they don't have mech coupling, do they?) then we have to sort out the coupling. Given my understanding of your physical problem, i can see 3 possible approaches:

1. it would be appropriate to simply set permeability=0 in the "masked" areas, without changing porosity at all in the MOOSE model (even if you know porosity does actually change). While you know that porosity may really change, what this method does is immobilize the water in those masked areas, and when/if they become unmasked then it can be ready to flow. While masked, temperature changes will cause porepressure changes in those masked regions - i'm not sure whether this will lead to poor convergence (eg, T increases might lead to enormous P increases, or T decreases might lead to P going negative) - you can calculate a few cases by hand as a guide (just keep density fixed in your analytical calculations, and see how P depends on T). When unmasked, you will suddenly create/destroy water if the porosity is not the same as it was upon being masked.

2. Probably a better alternative is to set porosity=0=permeability, which will result in MOOSE not calculating any fluid flow stuff in the masked areas. Your true porosity might be nonzero, but forcing porosity=0=permeability will mean MOOSE simply doesn't bother with any fluid stuff in the masked area (except for the fluid contributing to thermal mass and conductivity), so the masked areas won't contribute significantly to any compute time. This will keep P fixed, even if T changes. When unmasked, P will be what it was before masking, and you will suddenly create/destroy if the porosity is not the same as it was upon being masked. I think this is possibly the most fool-proof approach, but it depends on the computational method used to do the coupling (eg, with files, or something else).

3. An alternative, mentioned by @GiudGiud , is to simply have no "flow" Kernels on those masked areas. The problem with this is you have to define some mesh subdomains defining the masked areas, but it is the cleanest as far as MOOSE is concerned.

[garciapintado]

Hello Andy,
Thanks!
Yes, with SimpleFluidProperties, it is working fine in all tests.
Also, according to what you describe, I agree we'd then need a more clever extrapolation of the thermodynamic tables.
I did it first by extrapolating the slopes at the extremes, but by doing so the patterns were not consistently fine (as I felt about them), and opted for the 0-derivative approach.

For the coupling with our external [matlab] mechanical code, my plan is to pass the dynamic mechanical pressure [and plastic strain rate [only the IInd invariant], fields to MOOSE. Then add some [minor] code to Materials to include the parametrization I already have in Matlab in MOOSE. That is, it would mostly be based on PorousFlowPorosity and an additional option for permeability to be applied after Kozeny-Carman permeability-porosity. Scalar nominal parameters for porosity and permeability, would be specific for each subdomain. At the moment I neglect the feedback from the porepressure on the mechanical [i.e. the only feedback into the mechanical is thermal]. If the mechanical needs remeshing [which happens now and then, as the mesh deforms], the porepressure and thermal fields, which are passed and back and forth to the hydrothermal code, are just remeshed alongside the rest of mechanical variables [stress fields...]. This, in effect, creates/destroys water everytime that the hydrothermal part restarts, but as this happens every 10 kyr in our simulations, I think we can live with that for our big-scale purposes.

If we chose the option of 0 porosity 0 permeabilty for the "flow masked" areas [defined as special blocks], these blocks could just use a PorousFlowFlowPorosityConst, and PorousFlowPermeabilityConst instead in these zones.

Yes, the synthetic tests I sent have no mechanical coupling.

I'll do then some additional test with some variation in the extrapolation of the thermodynamic tables, and play with the point 2 you mention for coupling and will come back with that [if you have some algorithm at hand for extrapolating the EOS tables it would be great, but if not, I'll see to figure out something parsimonious but avoiding the zero-derivatives]

For the alternative 3 you describe, for which @GiudGiud sent a link to an example, I could automate the input [that is, creating some special blocks for these areas] and I feel computationally and physically is probably the best option. But it is obscure to me now how one would set up the boundary conditions. They would need to be along the boundaries of a subdomain for the advective part and it is currently confusing to me how one would organize the input file relating the individual kernels, blocks and BCs. I think we could not use the standard PorousFlowFulllySaturated action as it it, but rather rather substitute this action block in the input file, by a step-by-step input file, where it has to be made specific where each kernel applies to. Is this so?

I also need to do my homework on all this... allow me for a couple of days to become more familiar with the input file syntax.

[WilkAndy]

Sorry, i don't have any "extrapolator". I think i'd start with something simple. Eg, for density, r = r(P, T), i'd fit the existing data with r = a + bP + cT, then use that function to define the values at "infinity" (T=5000 or something), then use an inverse-distance weighted average for any intermediate points (eg, T=1000) if needed.

Before coding any C++ Materials, try passing porosity as a PiecewiseMultilinearFunction (or similar) defined by a file, then transferring to an AuxVariable, and passing that to a PorousFlowPorosityConst to define a constant porosity. I know that's not what you eventually want, but i feel if that works numerically and the file stuff works, then what you want will also work. Passing the new porosity in a file shouldn't produce any numerical problems, but i always try porosity changes first, since they are difficult numerically, because water is so extremely stiff, while permeability changes are much easier numerically.

If you decide to do alternative (3), then the MOOSE part won't be too difficult. You'll have to pass stuff from matlab to moose, so alternative (1) or (2) need to be done in some shape-or-form anyway, so probably good to get that sorted first. Aren't the BCs just no-flow (which translates to no BC stuff in your MOOSE input file)? It won't be a big deal unravelling PorousFlowFullySaturated into a whole lot of Kernels and Materials - i can help with that. But we'll have to be careful with recording porepressure in the "masked" areas. Your input file won't have any porepressure Variable living inside the masked areas, but when they become unmasked, you'll have to prescribe a value of porepressure there.

[garciapintado]

OK, I'll try first if with an updated, hopefully better, extrapolation of the EOS and fixed porosity the synthetic test for its most realistic scenario converges.

On alternative (3), I'm already passing stuff from matlab to MOOSE through EXODUS files [not in these synthetic but in the real case on which I initially failed to make it converge], prepared through an intermediate python script [called by matlab], using the package "exomerge" from Sandia NL, where I am passing full ICs files for temperature and porepressure, as well as nodesets for BCs ---I still haven't figured out how to prepare sidesets, but with nodeset BCs it is working for me by now---]. So, I believe all what could be needed as input could be automated as well.

On BCs my question is more about how MOOSE works when only some blocks are selected for a specific Kernel computation. That is, should the BCs be on nodes pertaining to these selected sub-blocks (either in their boundary or within their boundaries), or the BCs should be defined everywhere in the complete mesh [if so, in my mind I can just think off DirichletBC]. But don't bother in answering this [don't want to take an excess of your time], I'll do my homework with MOOSE. I would for sure need help on unravelling PorousFlowFullySaturated into the list of Kernels and Materials. But also, I'll start by studying it a bit more.

[WilkAndy]

If a Variable exists on some blocks only, then the remainder of the mesh is invisible to it. Eg, if P lives on some blocks then fluid will flow in those blocks only: as far as the fluid is concerned, it will see a domain with "holes" (the masked areas), so it will flow around those holes. The BCs should only be defined where the Variable exists: so if parts of your original boundary get masked, you'll have to create a new exodus boundary that consists of only the non-masked parts, and apply the DirichletBC there only.

I was pretty careful in describing all the Kernels and Materials added by the Action in the documentation https://mooseframework.inl.gov/source/actions/PorousFlowFullySaturated.html . Hopefully i didn't omit anything!

[garciapintado]

Hello Andy,
A "better" extrapolation of the thermodynamic tables has not reached a good end [I've tried several regression types, and some "parsimonious" extrapolation following the minimal slope in the borders of the original thermodynamic tables, and nothing seems to be able to solve the convergence problem when both high P and high T live together in the domain. Extrapolating with small slopes initially does not hit crazy values of temperature, but ultimately makes the run so slow that I had to stopped it. So, from everything I've tried still the 0-derivative extrapolation is the one which takes the simulations further.

So, I've gone for option [3] [mask the deepest block ---'mantle'---so no porepressure lives there.

I've now run the synthetic [4 layer on top of each other] indicating manually the Kernels and Material brought by the action PorousFlowFullySaturated [BTW, the documentation is also great!]. So, I've reproduced my former synthetic tests. Up to here, everything seems great.

Yet, now I'm stuck in how to specify that one variable [porepressure in this case] only lives in some blocks. I do not seem to find this in the documentation. So, I am getting a lot of errors, that I believe come because the porepressure Variable is seen everywhere:

Material property 'PorousFlow_saturation_qp', requested by 'lambda_mantle' is not defined on block mantle
Material property 'PorousFlow_saturation_qp', requested by 'lambda_mantle_neighbor' is not defined on block mantle
...
Material property 'dPorousFlow_porepressure_nodal_dvar', requested by 'PorousFlow_density_nodal_all' is not defined on block mantle
...
and so on. Is there any example where a variable only exists in some subDomain blocks? I really can't find it. Maybe I'm looking at the wrong places...

[WilkAndy]

By the way, who are you? You're obviously quite experienced with numerical modelling - it'd be nice to know where you're from, what you're researching, etc.

[garciapintado]

Thanks! Not terribly expert but trying to surf the waves :-)... I am Javier García-Pintado [Spanish], and work in MARUM (an excellence center in environmental marine sciences in Bremen, Germany, with a high percentage of resources going to experimental research [from biogeochemists to robotics people], and some resources for numerical analysis ---as my part---. I went into academia after some years in private industry in Spain (as mining engineer), and after a [more experimental than numerical] PhD in ephemeral rivers water quality, I moved to Reading (UK) to the Data Assimilation Research Center [DARC; most people there are mathematicians], and stayed there for 7 years working on assimilation of Synthetic Aperture Radar with the SWE [2D water shallow water equations] for flood forecast. However, for personal reasons [Brexit came to Europe, and then first my wife and then myself were offered to come to MARUM] we decided to move to Germany, and I had to change topics again. Now I work on the modelling of ocean basins within the MARUM cluster in collaboration with climate people, marine biogeochemists... The overarching goal in my part is to evaluate feedbacks between the formation of ocean basins at geological timescales [that's why we have our 2D viscoelastoplastic deformation model] and submarine hydrothermal systems [with implications for global reconstructions of heat release in relation natural climate variability, the initiation of life in cold-water submarine hydrothermal vents, black-smokers...]. So, after all these moves we're working in this great hybrid environment... AFAIK, what we are trying to do [now with your help] would become the first THM modelling at this broad scale for ocean basins.

My idea of trying MOOSE came because it clear that we need to go HPC [and MOOSE is so neatly done!]. Then, as geochemistry is added, the DOFs will highly increase so a proper scaling of the computing will be a must. I'm crossing fingers that this will work...

I've seen that in CSIRO, Mining Geomechanics, you are mostly mathematicians and mining engineers [which for me is a coincidence as I do not work formally as mining engineer, despite being one].
Cheers.

[garciapintado]

Hello Andy,
After looking a bit at the code and the options you mentioned I chose to attempt a simple modification to your own Kernel PorousFlowEnergyTimeDerivative as I saw that HeatCapacityConductionTimeDerivative does not include the density, and also under the hope that this would be more seamless integrated within PorousFlow (as temperature seems to be recognised a a PorousFlow variable?). Basically I added a boolean parameter neglect_porosity to the Kernel, which if false sets the number of phases to 0 and I saw that then you consider that there is not fluid present in the object.
So, I’m calling the version allowing this HeatCapacityConductionTimeDerivativeNopore for the mantle block with neglect_porosity=true.

As for PorousFlowHeatHeatConduction, I believe that now it does not generate any issue [I’ve just added a negligible PorousFlowPorosityConst porosity], and then a PorousFlowThermalConductivityFromPorosity for the mantle.

With this, plus making the Variable pore pressure available in all domain as well as initialising it in all the domain and adding the ‘mantle’ block to Material PorousFlow1PhaseFullySaturated, I’ve got to get rid of the initialisation errors.

But still, something I’ve done must be non-sensical, as I am getting the error:
“Linear solve did not converge due to DIVERGED_PC_FAILED iterations 0
PC failed due to FACTOR_OUTMEMORY”

Maybe the added neglect_porosity is not working as I expected?

Note that if I activate all blocks for all the kernels, the same simulation runs fine with the same preconditioning and Executioner.
I’m struggling to understand why this error, with all the environment [C++ & MOOSE] new to me…

Following the posting guidelines by @GiudGiud and guessing this is at this stage a convergence error, I am sending a link to GitHub repository, where I've put the files in case you could help: https://github.com/garciapintado/MOOSE_porous_flow_stuff
syn00.i : generates the mesh
syn08.i : all blocks active; converges very well
syn09.i : inactive Kernels for advective flow in the mantle: gives the error indicated above
PorousFlowEnergyTimeDerivativeNopore.h : exactly as yours but with extra parameter to neglect porosity
PorousFlowEnergyTimeDerivativeNopore.C: “ “ "

Have a nice weekend,
Javier

[garciapintado]

Hello Andy,
As in this setup the thermal advection kernel is only applied to some subdomain blocks (everywhere expect in the mantle), but the thermal conduction is applied everywhere, I am wondering if the issue may come from the thermal boundary conditions, which are specified for the top and bottom boundaries of the overall domain. Could it be that this is clashing with the subdomain application of the thermal advection?
Many thanks!
Javier

[WilkAndy]

The problem appears to be that mumps cannot invert Ax=b, which suggests your Jacobian is bad, or your mesh is bad. I'm not sure what's wrong. Can you try with the HeatCapacityConductionTimeDerivative Kernel instead? (I don't actually see the problem with your new Kernel, however.) Can you try to make the Variables order = first instead? It could be related to the TVD advection: can you try to turn that off (use full upwinding intead)?

[garciapintado]

Hello,
I’ve tried all the last indications [it would be very unlikely that the mesh has any issue as when all blocks are considered it works fine].
Note: I am not sure about how to add the input of HeatCapacityConductionTimeDerivative [as there is no example, and HeatCapacity is an extensive property] , but I’ve seen a SpecificHeatConductionTimeDerivate which I understand [and I think the intent of using this is what you mean?].

So, I’ve tried the SpecificHeatConductionTimeDerivative Kernel in the mantle, and making the variables LAGRANGE; 1st order as well. I’ve also tried the full-unwinding.
None of the modifications, one-by-one, seems to work. I am pasting here the input file, in case it is useful. As it is, including that the last modifications, the initialisation goes through, but the convergence error persists.

I’ve also done an alternative attempt, which I’ll comment in a separate post to avoid mixing details

[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
xmin = -200000
xmax = 200000
ymin = -150000
ymax = 0
elem_type = TRI3
nx = 100
ny = 75
[]
[lcrust]
type = SubdomainBoundingBoxGenerator
input = 'gen'
block_id = 1
#block_name = lcrust
bottom_left = '-200000 -40000 0'
top_right = '200000 -20000 0'
[]
[ucrust]
type = SubdomainBoundingBoxGenerator
input = 'lcrust'
block_id = 2
#block_name = ucrust
bottom_left = '-200000 -20000 0'
top_right = '200000 -10000 0'
[]
[unsat]
type = SubdomainBoundingBoxGenerator
input = 'ucrust'
block_id = 3
#block_name = unsat
bottom_left = '-200000 -10000 0'
top_right = '200000 0 0'
[]
[rename]
type = RenameBlockGenerator
input = 'unsat'
old_block = '0 1 2 3'
new_block = 'mantle lcrust ucrust unsat'
[]
[topnopore]
type = SideSetsBetweenSubdomainsGenerator
input = 'rename'
new_boundary = 'top_nopore'
primary_block = 'mantle'
paired_block = 'lcrust'
[]
[]

[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 -9.81 0' # make global, so no need to add it to the userobject requiring it
[]

[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure temperature'
number_fluid_phases = 1
number_fluid_components = 1
[]
[]

[Variables]
[porepressure]
block = 'mantle lcrust ucrust unsat' # blocks where this variable exists
family = LAGRANGE
order = FIRST
[]
[temperature]
family = LAGRANGE
order = FIRST
scaling = 1E-08
[]
[]

[Kernels] # as added by the action [PorousFlowFullySaturated]
[pp_time_derivative]
type = PorousFlowMassTimeDerivative
block = 'lcrust ucrust unsat'
variable = porepressure
[]
[pp_upwind_advectiveflux_kernel]
type = PorousFlowFullySaturatedAdvectiveFlux
block = 'lcrust ucrust unsat'
variable = porepressure
[]
[heat_time_derivative]
type = PorousFlowEnergyTimeDerivative
block = 'lcrust ucrust unsat'
variable = temperature
[]
[heat_time_derivative_nopore]
type = SpecificHeatConductionTimeDerivative
block = 'mantle'
variable = temperature
lumping = true
density = density_mantle
specific_heat = specific_heat_mantle
[]
[heat_conduction]
type = PorousFlowHeatConduction
variable = temperature
[]
[heat_upwind_advectiveflux_kernel]
type = PorousFlowFullySaturatedUpwindHeatAdvection
block = "lcrust ucrust unsat"
variable = temperature
[]
[]

[ICs]
[porepressure_IC]
type = FunctionIC
#block = "lcrust ucrust unsat"
variable = porepressure # [Pa]
function = '9.81*(-y)*1000'
[]
[temperature_IC]
type = FunctionIC
variable = temperature
function = '273.15+5+(-y)*600/150000' # => 878.15 = 273.15+5+600 at the bottom
[]
[]

[BCs]
[Ptop]
type = DirichletBC
variable = porepressure
value = 101325.0 # [Pa] 1 atm
boundary = top
[]
[Ttop]
type = DirichletBC
variable = temperature
value = 278.15 #
boundary = top
[]
[Tbot]
type = DirichletBC
variable = temperature
value = 878.15
boundary = top_nopore
[]
[]

[Modules]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
viscosity = 1.0E-3
density0 = 1000.0
[]
[]
[]

[Materials]
[materials_mantle]
type = GenericConstantMaterial
block = mantle
prop_names = 'thermal_conductivity_mantle specific_heat_mantle density_mantle'
prop_values = '3.3 1200.0 3360.0'
[]
[porosity_mantle]
type = PorousFlowPorosityConst # 'at_nodes=False' by default
block = mantle
porosity = 0.03
[]
[porosity_lcrust]
type = PorousFlowPorosity
block = lcrust
porosity_zero = 0.05
thermal = true
thermal_expansion_coeff = 3.66E-05
reference_temperature = 273.15
fluid = true
solid_bulk = 10.E09
biot_coefficient = 0.8
biot_coefficient_prime = 0.8
mechanical = false
[]
[porosity_ucrust]
type = PorousFlowPorosity
block = ucrust
porosity_zero = 0.05
thermal = true
thermal_expansion_coeff = 3.38E-05
reference_temperature = 273.15
fluid = true
solid_bulk = 10.E09
biot_coefficient = 0.8
biot_coefficient_prime = 0.8
mechanical = false
[]
[porosity_unsat]
type = PorousFlowPorosity
block = unsat
porosity_zero = 0.05
thermal = true
thermal_expansion_coeff = 3.38E-05
reference_temperature = 273.15
fluid = true
solid_bulk = 10.E09
biot_coefficient = 0.8
biot_coefficient_prime = 0.8
mechanical = false
[]
[permeability_mantle]
block = mantle
type = PorousFlowPermeabilityConst # 'at_nodes=False' by default
permeability = '1.0E-18 0 0 0 1.0E-18 0 0 0 1.0E-18'
[]
[permeability_lcrust]
type = PorousFlowPermeabilityKozenyCarman
block = lcrust
n = 3
m = 2
phi0 = 0.05
k0 = 9.0E-16
poroperm_function = 'kozeny_carman_phi0'
[]
[permeability_ucrust]
type = PorousFlowPermeabilityKozenyCarman
block = ucrust
n = 3
m = 2
phi0 = 0.05
k0 = 1.0E-15
poroperm_function = 'kozeny_carman_phi0'
[]
[permeability_unsat]
type = PorousFlowPermeabilityKozenyCarman
block = unsat
n = 3
m = 2
phi0 = 0.05
k0 = 1.0E-22
poroperm_function = 'kozeny_carman_phi0'
[]
[internal_energy_mantle] # 'at_nodes=True' by default
type = PorousFlowMatrixInternalEnergy # this Material calculated the internal energy of solid rock grains
block = mantle # internal_energy [J.K-1.m-3] = density*specific_heat_capacity
density = 3360.0 # [kg.m-3] density of rock grains
specific_heat_capacity = 1200.0 # [J.kg-1.K-1] specific heat capacity of rock grains
[]
[internal_energy_lcrust]
type = PorousFlowMatrixInternalEnergy # this Material calculated the internal energy of solid rock grains
block = lcrust
density = 2850.0 # [kg.m-3] density of rock grains
specific_heat_capacity = 1200.0 # [J.kg-1.K-1] specific heat capacity of rock grains
[]
[internal_energy_ucrust]
type = PorousFlowMatrixInternalEnergy # this Material calculated the internal energy of solid rock grains
block = ucrust
density = 2700.0 # [kg.m-3] density of rock grains
specific_heat_capacity = 1200.0 # [J.kg-1.K-1] specific heat capacity of rock grains
[]
[internal_energy_unsat]
type = PorousFlowMatrixInternalEnergy # this Material calculated the internal energy of solid rock grains
block = unsat
density = 2700.0 # [kg.m-3] density of rock grains
specific_heat_capacity = 1200.0 # [J.kg-1.K-1] specific heat capacity of rock grains
[]
[lambda_mantle]
type = PorousFlowThermalConductivityFromPorosity # rock-fluid combined thermal conductivity by weighted sum of rock and fluid conductivities
block = mantle
lambda_f = '3.3 0 0 0 3.3 0 0 0 3.3'
lambda_s = '3.3 0 0 0 3.3 0 0 0 3.3'
[]
[lambda_lcrust]
type = PorousFlowThermalConductivityFromPorosity
block = lcrust
lambda_f = '0.56 0 0 0 0.56 0 0 0 0.56'
lambda_s = '2.5 0 0 0 2.5 0 0 0 2.5'
[]
[lambda_ucrust]
type = PorousFlowThermalConductivityFromPorosity
block = ucrust
lambda_f = '0.56 0 0 0 0.56 0 0 0 0.56'
lambda_s = '2.3 0 0 0 2.3 0 0 0 2.3'
[]
[lambda_unsat]
type = PorousFlowThermalConductivityFromPorosity
block = unsat
lambda_f = '0.56 0 0 0 0.56 0 0 0 0.56'
lambda_s = '2.3 0 0 0 2.3 0 0 0 2.3'
[]

[porepressure_material]
type = PorousFlow1PhaseFullySaturated # at_nodes=false by default
block = 'mantle lcrust ucrust unsat'
porepressure = porepressure
[]
[temperature_material]
type = PorousFlowTemperature # at_nodes=false by default
temperature = temperature
[]
[massfrac]
type = PorousFlowMassFraction # at_nodes=false by default
block = 'lcrust ucrust unsat' # list of blocks where this object applies to
# mass_fraction_vars # nod needed when num_phases=num_component
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid # see documentation for this material regarding the choice of units
block = 'lcrust ucrust unsat'
fp = the_simple_fluid # this Material is at_nodes=false by default
phase = 0
[]
[simple_fluid_noflow]
type = PorousFlowSingleComponentFluid # see documentation for this material regarding the choice of units
block = 'mantle'
fp = the_simple_fluid # this Material is at_nodes=false by default
phase = 0
[]
[effective_fluid_pressure] # create effective fluid pressure [requested by PorousFlowPorosity even it has not mechanical coupling]
block = 'lcrust ucrust unsat'
type = PorousFlowEffectiveFluidPressure
[]
[nearest_qp]
type = PorousFlowNearestQp # atNodes=false by default
block = 'mantle lcrust ucrust unsat' # should mantle be neglected here?
[]
[relperm] # required by PorousFlowDarcyVelocityComponent AuxKernels
type = PorousFlowRelativePermeabilityConst # atNodes=false by default
block = 'mantle lcrust ucrust unsat' # ERROR when mantle not specified
phase = 0
kr = 1
[]
[]

[Preconditioning]
active = smp_lu_mumps
[smp_lu_mumps]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu mumps'
[]
[]

[Executioner]
type = Transient
solve_type = Newton
end_time = 315576000000 # [s] 10000 year
#end_time = 1E6 # [s] ~[12 days]
dtmax = 3.2E10 # [s] ~1000 year. advanced parameter. Maximum timestep size in an adaptive run. Default to 1E30
nl_max_its = 25 # solver parameter. Max Nonlinear Iterations. Default to 50
l_max_its = 100 # solver parameter. Max Linear Iterations. Default to 10000
nl_abs_tol = 1E-06 # solver parameter. Nonlinear absolute tolerance. Default to 1E-50
nl_rel_tol = 1E-08 # solver parameter. Nonlinear Relative Tolerance. Default to 1E-08
scheme = 'implicit-euler' # the default TimeIntegrator
[TimeStepper] # TimeStepper subsystem [block always nested within the Executioner block]
type = IterationAdaptiveDT # adjust the timestep based on the number of iterations
optimal_iterations = 10 # target number of nonlinear iterations
dt = 1E5 # ~[27 h]
growth_factor = 2
cutback_factor = 0.5
[]
[]

[Outputs]
[ou]
type = Exodus
#output_material_properties = true
[]
[]

[garciapintado]

Hello,
And this is an alternative attempt. Here I thought about constraining the variables controlled by the PorousFlow dictator {porepressure, temperature} to the PorousFlow subdomains (everywhere but the mantle). That is, neither in their definition nor everywhere else they are supposed to exist in the mantle. None of the related Kernels or Materials makes any reference to the mantle [at least not on purpose].

Then I’ve created and alternative variable, non-controlled by the PorousFlow dictator, named {temperature_mantle} which only exists in the mantle, and is subject only to conductive evolution, through the use of non PorousFlow Kernels: SpecificHeatConductionTimeDerivative and HeatConduction [after your previous indications], with their specific non-PorousFlow parameters.
Also, there is now an additional internal side between the mantle and the above domain (the one that house be visible by the PorousFlow dictator). And there are boundary conditions specific to temperature and to temperature_mantle. Step-by-step, to remove any possibility of creating an error because of any interaction between temperature and temperature_mantle, they are not coupled. I would not know, anyway, how to coupled them yet, but I guess this should be possible somehow through their common boundary. So, I’d just expect them to evolve as independent variables.

That is: the aim of the test is [a] to see if PorousFlow can work only on a subdomain(s), and [b] if this worked, then I'd expect the overall “multiphysics” problem could be approached with the current PorousFlow+MOOSE [i.e. only heat conduction in the mantle, and the THC —thermo-hydro-chemical problem everywhere else—], provided that later on temperatureand temperature_mantle can be coupled, so their coupling would substitute the values at their common boundary.

Well, unfortunately, PorousFlow gives an error again at initialisation of Materials:

…
Material property 'dPorousFlow_porepressure_nodal_dvar', requested by 'PorousFlow_density_nodal_all' is not defined on block mantle
Material property 'dPorousFlow_saturation_nodal_dvar', requested by 'PorousFlow_density_nodal_all' is not defined on block mantle
Material property 'dPorousFlow_temperature_nodal_dvar', requested by 'PorousFlow_density_nodal_all' is not defined on block mantle
…

It appears that at initialisation, PorousFlowDictator is forcing some Action to occur everywhere in the domain [perhaps PorousFlowAddMaterialAction or PorousFlowAddMaterialJoiner?]. If this is so, it seems that as it stands, it can’t be constrained to work in specific subdomains, preventing it from multiphysics coupling with adjacent domains. Could this be so?

If so, and there is anything I could do to help adding this possibility [facing my severe lack of C++ knowledge], please let me know.

I am pasting the self-contained input file for this test.

[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
xmin = -200000
xmax = 200000
ymin = -150000
ymax = 0
elem_type = TRI6
nx = 100
ny = 75
[]
[lcrust]
type = SubdomainBoundingBoxGenerator
input = 'gen'
block_id = 1
bottom_left = '-200000 -40000 0'
top_right = '200000 -20000 0'
[]
[ucrust]
type = SubdomainBoundingBoxGenerator
input = 'lcrust'
block_id = 2
bottom_left = '-200000 -20000 0'
top_right = '200000 -10000 0'
[]
[unsat]
type = SubdomainBoundingBoxGenerator
input = 'ucrust'
block_id = 3
bottom_left = '-200000 -10000 0'
top_right = '200000 0 0'
[]
[rename]
type = RenameBlockGenerator
input = 'unsat'
old_block = '0 1 2 3'
new_block = 'mantle lcrust ucrust unsat'
[]
[topnopore]
type = SideSetsBetweenSubdomainsGenerator
input = 'rename'
new_boundary = 'top_nopore'
primary_block = 'mantle'
paired_block = 'lcrust'
[]
[]

[GlobalParams]
PorousFlowDictator = dictator
gravity = '0 -9.81 0' # make global, so no need to add it to the userobject requiring it
[]

[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure temperature'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pp_KT_advectiveflux_onecomp_userobj]
type = PorousFlowAdvectiveFluxCalculatorSaturated # [1-phase, 1-comp, fully saturated] one kernel for each fluid component
block = 'lcrust ucrust unsat' # OPTIONAL, list of blocks that this object will be applied to
flux_limiter_type = superbee
multiply_by_density = true # default. Implies the advective flux is multiplied by density, so it is a mass flux
[]
[heat_KT_advectiveflux_userobj]
type = PorousFlowAdvectiveFluxCalculatorSaturatedHeat
block = 'lcrust ucrust unsat'
flux_limiter_type = superbee
multiply_by_density = true # default. Implies the advective flux is multiplied by density, so it is a mass flux
[]
[]

[Variables]
[porepressure]
block = 'lcrust ucrust unsat' # blocks where this variable exists
family = LAGRANGE
order = SECOND
[]
[temperature] # [K] parsed initial condition in the ICs block below
block = 'lcrust ucrust unsat' # blocks where this variable exists
family = LAGRANGE
order = SECOND
scaling = 1E-08 # this variable scaling brings the residual R_T to the same order of magnitude that the Pressure residual R_P
[]
[temperature_mantle] # [K] parsed initial condition in the ICs block below
block = 'mantle' # blocks where this variable exists
family = LAGRANGE
order = SECOND
scaling = 1E-08 # this variable scaling brings the residual R_T to the same order of magnitude that the Pressure residual R_P
[]
[]

[Kernels] # as added by the action [PorousFlowFullySaturated]
[pp_time_derivative] # one kernel for each fluid component
type = PorousFlowMassTimeDerivative # these kernels lump the fluid-component mass to the nodes to ensure superior numerical stabilization
block = 'lcrust ucrust unsat'
variable = porepressure
[]
[pp_KT_advectiveflux_kernel]
type = PorousFlowFluxLimitedTVDAdvection # one kernel for each fluid component
block = 'lcrust ucrust unsat'
advective_flux_calculator = pp_KT_advectiveflux_onecomp_userobj # PorousFlowAdvectiveFluxCalculator UserObjectName
variable = porepressure
[]
[heat_time_derivative]
type = PorousFlowEnergyTimeDerivative # this kernel lumps the heat energy-density to the nodes to ensure superior numerical stabilization
block = 'lcrust ucrust unsat'
variable = temperature
# save_in = # name of auxiliary variable to save this Kernel residual contributions to.
[]
[heat_conduction]
type = PorousFlowHeatConduction
block = 'lcrust ucrust unsat'
variable = temperature
[]
[heat_KT_advectiveflux_kernel]
type = PorousFlowFluxLimitedTVDAdvection # one kernel for each fluid component
block = "lcrust ucrust unsat"
advective_flux_calculator = heat_KT_advectiveflux_userobj # PorousFlowAdvectiveFluxCalculator UserObjectName
variable = temperature
[]
[heat_time_derivative_mantle]
type = SpecificHeatConductionTimeDerivative
block = 'mantle'
variable = temperature_mantle
lumping = true
density = density_mantle
specific_heat = specific_heat_mantle
[]
[heat_conduction_mantle]
type = HeatConduction
block = 'mantle'
variable = temperature_mantle
diffusion_coefficient = thermal_conductivity_mantle
[]
[]

[ICs]
[porepressure_IC]
type = FunctionIC
block = "lcrust ucrust unsat"
variable = porepressure # [Pa]
function = '9.81*(-y)*1000'
[]
[temperature_IC]
type = FunctionIC
block = "lcrust ucrust unsat"
variable = temperature
function = '273.15+5+(-y)600/150000' # => 878.15 = 273.15+5+600 at the bottom [273.15+5+600/15000040000]
[]
[temperature_mantle_IC]
type = FunctionIC
block = "mantle"
variable = temperature_mantle
function = '273.15+5+(-y)600/150000' # => 878.15 = 273.15+5+600 at the bottom [273.15+5+600/15000040000]
[]
[]

[BCs]
[Ptop]
type = DirichletBC # this is a NodalBC
variable = porepressure
value = 101325.0 # [Pa] 1 atm
boundary = top
[]
[Ttop]
type = DirichletBC
variable = temperature
value = 278.15 #
boundary = top
[]
[Tbot]
type = DirichletBC
variable = temperature
value = 438.15 # matching grad_T_y = 600/150000
boundary = top_nopore
[]
[Ttop_mantle]
type = DirichletBC
variable = temperature_mantle
value = 438.15 # matching grad_T_y = 600/150000
boundary = top_nopore
[]
[Tbot_mantle]
type = DirichletBC
variable = temperature_mantle
value = 878.15 # matching grad_T_y = 600/150000
boundary = bottom
[]
[]

[Modules]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
viscosity = 1.0E-3
density0 = 1000.0
[]
[]
[]

[Materials]
[materials_mantle]
type = GenericConstantMaterial
block = mantle
prop_names = 'thermal_conductivity_mantle specific_heat_mantle density_mantle'
prop_values = '3.3 1200.0 3360.0'
[]
[porosity_lcrust]
type = PorousFlowPorosity
block = lcrust
porosity_zero = 0.05
thermal = true
thermal_expansion_coeff = 3.66E-05
reference_temperature = 273.15
fluid = true
solid_bulk = 10.E09
biot_coefficient = 0.8
biot_coefficient_prime = 0.8
mechanical = false
[]
[porosity_ucrust]
type = PorousFlowPorosity
block = ucrust
porosity_zero = 0.05
thermal = true
thermal_expansion_coeff = 3.38E-05
reference_temperature = 273.15
fluid = true
solid_bulk = 10.E09
biot_coefficient = 0.8
biot_coefficient_prime = 0.8
mechanical = false
[]
[porosity_unsat]
type = PorousFlowPorosity
block = unsat
porosity_zero = 0.05
thermal = true
thermal_expansion_coeff = 3.38E-05
reference_temperature = 273.15
fluid = true
solid_bulk = 10.E09
biot_coefficient = 0.8
biot_coefficient_prime = 0.8
mechanical = false
[]
[permeability_lcrust]
type = PorousFlowPermeabilityKozenyCarman
block = lcrust
n = 3
m = 2
phi0 = 0.05
k0 = 9.0E-16
poroperm_function = 'kozeny_carman_phi0'
[]
[permeability_ucrust]
type = PorousFlowPermeabilityKozenyCarman
block = ucrust
n = 3
m = 2
phi0 = 0.05
k0 = 1.0E-15
poroperm_function = 'kozeny_carman_phi0'
[]
[permeability_unsat]
type = PorousFlowPermeabilityKozenyCarman
block = unsat
n = 3
m = 2
phi0 = 0.05
k0 = 1.0E-22
poroperm_function = 'kozeny_carman_phi0'
[]
[internal_energy_lcrust]
type = PorousFlowMatrixInternalEnergy
block = lcrust
density = 2850.0 # [kg.m-3] density of rock grains
specific_heat_capacity = 1200.0 # [J.kg-1.K-1] specific heat capacity of rock grains
[]
[internal_energy_ucrust]
type = PorousFlowMatrixInternalEnergy
density = 2700.0
specific_heat_capacity = 1200.0
[]
[internal_energy_unsat]
type = PorousFlowMatrixInternalEnergy
block = unsat
density = 2700.0
specific_heat_capacity = 1200.0
[]
[lambda_lcrust]
type = PorousFlowThermalConductivityFromPorosity
block = lcrust
lambda_f = '0.56 0 0 0 0.56 0 0 0 0.56'
lambda_s = '2.5 0 0 0 2.5 0 0 0 2.5'
[]
[lambda_ucrust]
type = PorousFlowThermalConductivityFromPorosity
block = ucrust
lambda_f = '0.56 0 0 0 0.56 0 0 0 0.56'
lambda_s = '2.3 0 0 0 2.3 0 0 0 2.3'
[]
[lambda_unsat]
type = PorousFlowThermalConductivityFromPorosity # rock-fluid combined thermal conductivity by weighted sum of rock and fluid conductivities
block = unsat
lambda_f = '0.56 0 0 0 0.56 0 0 0 0.56'
lambda_s = '2.3 0 0 0 2.3 0 0 0 2.3'
[]
[porepressure_material]
type = PorousFlow1PhaseFullySaturated # at_nodes=false by default
block = 'lcrust ucrust unsat'
porepressure = porepressure
[]
[temperature_material]
type = PorousFlowTemperature # at_nodes=false by default
block = 'lcrust ucrust unsat' # I think I would not need this here, as temperature exists in all blocks
temperature = temperature
[]
[massfrac]
type = PorousFlowMassFraction # at_nodes=false by default
block = 'lcrust ucrust unsat' # list of blocks where this object applies to
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid # see documentation for this material regarding the choice of units
block = 'lcrust ucrust unsat'
fp = the_simple_fluid # this Material is at_nodes=false by default
phase = 0
[]
[effective_fluid_pressure] # create effective fluid pressure [is requested by PorousFlowPorosity even it has not mechanical coupling]
block = 'lcrust ucrust unsat'
type = PorousFlowEffectiveFluidPressure
[]
[nearest_qp]
type = PorousFlowNearestQp # atNodes=false by default
block = 'lcrust ucrust unsat' # should mantle be neglected here?
[]
[relperm] # required by PorousFlowDarcyVelocityComponent AuxKernels
type = PorousFlowRelativePermeabilityConst # atNodes=false by default
block = 'lcrust ucrust unsat' # ERROR when mantle not specified
phase = 0
kr = 1 # default, anyway
[]
[]

[Preconditioning]
active = smp_lu_mumps
[smp_lu_mumps]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu mumps'
[]
[]

[Executioner]
type = Transient
solve_type = Newton
end_time = 315576000000 # [s] 10000 year
#end_time = 1E6 # [s] ~[12 days]
dtmax = 3.2E10 # [s] ~1000 year. advanced parameter. Maximum timestep size in an adaptive run. Default to 1E30
nl_max_its = 25 # solver parameter. Max Nonlinear Iterations. Default to 50
l_max_its = 100 # solver parameter. Max Linear Iterations. Default to 10000
nl_abs_tol = 1E-06 # solver parameter. Nonlinear absolute tolerance. Default to 1E-50
nl_rel_tol = 1E-08 # solver parameter. Nonlinear Relative Tolerance. Default to 1E-08
scheme = 'implicit-euler' # the default TimeIntegrator
[TimeStepper] # TimeStepper subsystem [block always nested within the Executioner block]
type = IterationAdaptiveDT # adjust the timestep based on the number of iterations
optimal_iterations = 10 # target number of nonlinear iterations
dt = 1E5 # ~[27 h]
growth_factor = 2
cutback_factor = 0.5
[]
[]

[Outputs]
[ou]
type = Exodus
#output_material_properties = true
[]
[]

[WilkAndy]

Hi @garciapintado . This looks like a bug to me. We need to look into it.

[WilkAndy]

Hi @cpgr , i think we need to fix a bug in PorousFlow . The problem is when the PorousFlow Variables are block restricted. MOOSE then complains it needs certain Materials on blocks where the PorousFlow Variables do not exist. For instance, if the PorousFlow Variables (and Materials, etc) do not appear on block = mantle, then MOOSE spits back things like:

Material property 'dPorousFlow_porepressure_nodal_dvar', requested by 'PorousFlow_density_nodal_all' is not defined on block mantle
Material property 'dPorousFlow_porepressure_nodal_dvar', requested by 'PorousFlow_density_nodal_all_face' is not defined on block mantle
Material property 'dPorousFlow_porepressure_nodal_dvar', requested by 'PorousFlow_density_nodal_all_neighbor' is not defined on block mantle

Note the _face and _neighbor in the last two messages. What is this?!

(It appears that the "core" variables, such as PorousFlow_porepressure_nodal are not needed, but the derivatives dPorousFlow_porepressure_nodal_dvar are, but there are quite a few error messages, so i could just be missing something.)

Hopefully today i'll have time to create a simple example.

[garciapintado]

Hi @WilkAndy & @cpgr, Many thanks for looking into this!

[WilkAndy]

Hi @cpgr . I believe the problem is in the Joiners. As far as i can see PorousFlowAddMaterialJoiner.C does not take into account block-restricted materials. I think every time we addJoiner we'll have to pass block-restrictions along. Do you agree with this? I'm going to start coding it now, but would be good to have your confirmation.

[WilkAndy]

@garciapintado , i've got a fix for this. It will take some time for it to appear in MOOSE: perhaps 1 week? If you are in a hurry i can tell you how to get my code.

[garciapintado]

How fast! Many thanks! I think I can wait for a week to just do the git pull and in the meantime try to learn how to add an extra material myself through an Exodus input file [a REAL volumetric strain field which would replace the one currently taken from the mechanical coupling for PorousFlowPorosity, plus another REAL field to be used as factor for the permeability]. Basic stuff, but my first attempt to do some actual C++ MOOSE coding, and I have to do it anyway.

Shall we close this discussion as the masking topic is answered, or keep it open to follow up further steps in this MOOSE+PorousFlow application?

[WilkAndy]

Thanks for finding the bug, @garciapintado !

[cpgr]

Oh, I wasn't expecting this!

So if I understand this, the joiner makes a material everywhere, even where a variable and subsequent nodal material aren't even present? It all works when there is the same material defined on different blocks, as long as it is on all blocks, but fails when one block doesn't have that material (like in the example above).

Hopefully you can propagate the block info easily enough.

[WilkAndy]

Yea. Maybe i'm being lazy and not querying all the Materials multiple times. What do you think?

[WilkAndy]

Anyway, i implemented the Dictator idea and it works, so i know the Joiners were the problem

[cpgr]

This is probably the easiest fix.

[WilkAndy]

Heheheh, now i'm thinking i'm being lazy!!

[cpgr]

Efficient sounds better

[WilkAndy]

Fixing this problem in #20260

[WilkAndy]

Yea, good, let's close

[WilkAndy]

Hi @garciapintado - how is your modelling going?

[garciapintado]

Hi @WilkAndy, I believe it should work for me now!
A simple 2D test works as expected [temperature only controlled by diffusion in the mantle, and Porous flow for the top layers, and porepressure only existing in the top layers]. A bit more complicated tests are crashing but this must be my mistake. So I'll try to figure out the errors and hopefully come back with something more realistic running...

[garciapintado]

Hi @WilkAndy. Just to say that I've got a "realistic" synthetic domain running with tabulated fluid properties [water97]. Now, as the [synthetic] 'mantle' is not considered, the high pressure and temperature values that led to the previous crashes are not there anymore. This is great news to me!
I'll go on now with the real case simulations. One note is that I am using constant porosity and permeability [haven't managed to make it run with PorousFlowPorosity]. But in the application I'll pass constant fields of porosity & permeability pre-calculated in Matlab as 1st step to circumvent this...

[garciapintado]

Hi @WilkAndy. Only to report that I later noticed that geothermal gradients were not realistic in my synthetic tests. For geothermal gradients going to about 800 ºC at 120km depth the runs worked well, but when I changed that to a more standard [~1300ºC at 120 km depth]. The solver did not converge at all, often with warnings about water temperature being below the minimum for the EOS [eventually I got rid of these warning by extrapolating ---with 0-slope the tabulated EOS also towards lower temperature values, but still this did not solve the convergence problems]. The limiting gradient for converging runs was that to reach ~900ºC at 120 km depth].

So, I've been struggling with that [plus another matlab settings for the coupling]. In the end I've just now tried to switch from the Kuzmin-Turek stabilization to a full upwinding for both fluid pressure and temperature, and it turns out that now the runs converge very well for the more realistic geothermal gradient.

So, this is just to let you know that the full upwinding schemes converge better than KZ in these synthetic tests. I'll go now, finally, for the real case scenario with input from matlab... I'd presume that the more complicated new part for me will be to set the BCs, as some areas are under the ocean [for which Dirichlet are not OK for outflow pressures and temperatures --- the venting areas---], and some other are actually subaerial in the sides of the domain [actually I am still unsure if it is possible to go ahead with fully saturated flow, which I'd like if possible]. Anyway, I'll report results...

[garciapintado]

Sorry, just to note thatt KZ stabilisation with flux limiter VanLeer also converges very well (3 times slower than full upwinding but very fast still). So, the issue seemed to be the "superbee" flux limiter I was using before.

[WilkAndy]

OK, that's good news about the flux limiter, but .... i don't understand it! I'll have to have a think :-)