Use of "region" and "state" in PDE models with variable coefficients.

"state" and "region" sometimes show up in questions about solving PDEs with variable coefficients, but these things are basically undocumented. I am trying to model a 3-region heat-transfer system in which the thermal conductivity of one layer is a function of temperature (see https://www.mathworks.com/matlabcentral/answers/498820-how-to-write-anonymous-function-for-variable-coefficients-in-heat-transfer-problem). What I wrote in the other post does not give an error but also does not assign the proper value to the layer in question.

If anyone has a grasp of how to use state and region to assign variable properties in PDEs, as in the earlier post, please reply.

Answer 1

[Ravi Kumar]

Hi Allen,

Solver passes two input arguments to the function that you define, "region" and "state" are just place holder names. You can call them anything you want. Just be aware that the first argument, region, contains spatial data which you can use to compute k and second contains solution data to serve you the same purpose.

Now to your specific question on kFunc I think you are using the logical expression state.u<cractT to modify K. The input data, that solver sends to your function, state.u contains solution at several points withiin Face 3. Depending on your initial condition you start out with state.u<cractT yielding logical vector (logical zeros and ones). I am guessing this is not what you want. Also, it is hard to say if the solution actually crosses 900 to change the thermal conductivity. If the condition is never false, then you will always get k +k*(Nu-1) as thermal conductivity. I would suggest writing a full function, instead of anonymous function as:

function kOut = kFunc(region,state)
if any(isnan(state.u))
    kOut = nan(size(state.u));
    return;
end
kOut = k;
if any(state.u < crackT)
    kOut = k+k*(Nu-1)*ones(size(state.u));
end
end

Then specify, this kFunc's handle on to faces 3 as:

thermalProperties(thermalModel,'Face',3,'ThermalConductivity',@kFunc,'MassDensity',2500,'SpecificHeat',1000);

Regards,

Ravi

[Allen]

That is very instructive, thanks. I placed the function at the end of the script and replaced the thermalProperties line as specified above. Now the script fails with a "Function specifying a material property must accept 2 input arguments and return 1 output argument". I've run into this error at various times in this process. "region" is not specifically used in the function, but it is an argument.

[Allen]

More puzzling information about this issue. First, clarification. My goal is to set parts of Face 3 where T < crackT to k*Nu, and leave the rest at k (in other words, thermal conductivity goes up where the material cracks). Initial T of Face 3 is 1200. After much experimenting I have gone back to this code:

Nu = 10;
crackT = 1200;
kF = @(region,state) k+(k*(Nu-1)*double(state.u<crackT));

thermalProperties(thermalModel,'Face',1,'ThermalConductivity',k*Nu,'MassDensity',2500,'SpecificHeat',1000);
thermalProperties(thermalModel,'Face',2,'ThermalConductivity',k,'MassDensity',2500,'SpecificHeat',1000);
thermalProperties(thermalModel,'Face',3,'ThermalConductivity',kF,'MassDensity',2500,'SpecificHeat',1000);

My intent is that the kF function will set those nodes where T < crackT to k*Nu. This code runs without error. However, I have found that if crackT > 1200 then all of Face 3 is set to k*Nu, but if crackT <= 1200 then all of Face 3 remains at k. Thus, the reassignment is all-or-nothing and the logical test does nothing.

I have also tried the following function (modified from comment above) with "@kFunc" in the last line of code above, but get the "Function specifying a material property must accept 2...etc." error.

function kOut = kFunc(region,state)
if any(isnan(state.u))
    kOut = nan(size(state.u));
    return;
end
if state.u < crackT
    kOut = k*Nu*state.u;
else
    kOut = k;
end

Perhaps "region" should be specified somehow?

[Ravi Kumar]

I don't think you need to use region. But you do need to give the consideration that you have to send a vector of kOut that corresponds to each spactial point defined by (region.x,region.y). My be this is what you are tyring to do:

function kOut = kFunc(region,state)
if any(isnan(state.u))
    kOut = nan(size(state.u));
    return;
end
% Have a constatnt k vector for all points that solver is querying
kOut = k.*ones(size(region.x)); % You can use state.u to also here. 

if any(state.u < crackT)
    kOut(stat.u<crackT) = k*Nu*state.u; % Logical indexing into to points where you need different k
end
end

If this doesn's solve, can you provde the complete reproduction code, could be simplifed version of your problem? Then someone can reproduce ths issue and investigate further.

Regards,

Ravi

[Allen]

Here is a stripped-down version that gives the usual error:

Error using pde.ThermalMaterialAssignment/checkCoefFcnHdlArgCounts (line 368)

Function specifying a material property must accept 2

input arguments and return 1 output argument.

% Cooling of a planar sheet
% Define  variables

times = logspace(2,6,10);
secPerYear = 3.17e7;
dTdz = 0.02; % geothermal gradient in °C/m 
width = 10000; % width of model in m
mleft = 0; % left margin of intrusion
mright = 10000; % right margin of intrusion
depth = 30000; % depth of model in m
depth1 = -3000; % bottom of top layer
depth2 = -8000; % bottom of middle layer
depth3 = -30000; % bottom of lower layer

% Define rectangles and plot them to show edge designations
R1 = [3 4  0 width width 0 0  0 depth2 depth2]';
R2 = [3 4 mleft mright mright mleft depth2 depth2 depth1 depth1]';
R3 = [3 4  0 width width 0 depth3 depth3 depth2 depth2]';
gd = [R1 R2 R3];
sf = 'R1+R2+R3'; ns = [82 82 82;49 50 51]; g = decsg(gd,sf,ns); % I have no idea what this does
thermalModel = createpde('thermal','transient');
geometryFromEdges(thermalModel,g);
subplot(1,3,1)
pdegplot(thermalModel,'EdgeLabels','on','FaceLabels','on'); axis equal;

% Set boundary conditions
thermalBC(thermalModel,'Edge',1,'Temperature',0);
thermalBC(thermalModel,'Edge',3,'HeatFlux',0.04.*secPerYear);

% Generate and plot mesh
msh = generateMesh(thermalModel,'Hmax',1000);
subplot(1,3,2)
pdeplot(thermalModel);

% Set up properties of differential equation
% I am using time units of years (in thermal conductivity and time).
k = 2*secPerYear; % thermal conductivity in J/(a m K)
rho = 2500; % density in kg/m3

Nu = 1;
crackT = 900;

thermalProperties(thermalModel,'Face',1,'ThermalConductivity',k*Nu,'MassDensity',2500,'SpecificHeat',1000);
thermalProperties(thermalModel,'Face',2,'ThermalConductivity',k,'MassDensity',2500,'SpecificHeat',1000);
thermalProperties(thermalModel,'Face',3,'ThermalConductivity',@kFunc,'MassDensity',2500,'SpecificHeat',1000);

T0f1 = @(location)(-dTdz/Nu)*location.y; % this sets the initial conditions for face 1
thermalIC(thermalModel,T0f1,'Face',1);
thermalIC(thermalModel,1200,'Face',3); %overwrite for magma body
T0f2 = @(location)(-dTdz/Nu)*depth2-dTdz*(location.y-depth2); % this sets the initial conditions for face 2
thermalIC(thermalModel,T0f2,'Face',2);

% Solve the equations
times = logspace(0,6,50);
thermalresults=solve(thermalModel,times); %

% Plot results
T = thermalresults.Temperature; % extract temperature

subplot(1,3,3)
for i = 1:numel(times)
    pdeplot(thermalModel,'XYData',T(:,i),'Contour','on','Levels',0:100:1000)
    colormap(parula(200))
    xlim([0 width])
    ylim([depth3 0]);
    caxis([0 1000]);
    drawnow
end

% Auxiliary function for thermal conductivity
function kOut = kFunc(region,state)
if any(isnan(state.u))
    kOut = nan(size(state.u));
    return;
end
% Have a constant k vector for all points that solver is querying
kOut = k.*ones(size(region.x)); % You can use state.u to also here. 
if any(state.u < crackT)
    kOut(state.u<crackT) = k*Nu*state.u; % Logical indexing into to points where you need different k
end
end

[Ravi Kumar]

Thanks for the reproduction code. There are a couple mistake in kFunc. I fixed those in the following version:

function kOut = kFunc(region,state)
if any(isnan(state.u))
    kOut = nan(size(state.u));
    return;
end
% Have a constant k vector for all points that solver is querying
k = 2*3.17e7;
crackT = 900;
Nu = 1;
kOut = k.*ones(size(region.x)); % You can use state.u to also here. 
if any(state.u < crackT)
    kOut(state.u<crackT) = k*Nu*state.u(state.u<crackT); % Logical indexing into to points where you need different k
end
end

If you use this function, you won't get any setup related error. However, the solution does not converge. I get the warning:

Failure at t=2.536495e+03. Unable to meet integration tolerances without reducing

the step size below the smallest value allowed (7.275958e-12) at time t.

This means time-integrator did not converge at the notes time-step. This can happen if either solution does not exisit (nonlinear problem) or there is a some problems in the parameters provided as input.

Regards.

Ravi

[Allen]

Ravi,

Thanks very much for working on this. A few things:

1. Last line of code should be kOut(state.u<crackT) = Nu*kOut(state.u<crackT);
2. For some reason if Nu is not defined in kFunc, then the "2 input arguments" error happens.
3. There are some interesting interactions between Nu and crackT that cause the convergence error. I am experimenting with that and playing with tolerances. Sometimes it does not fail but just sits on the solve step for a long time.
4. I still do not understand the "2 input arguments" error, but at least things seem to be working! Thanks again. I hope this thread helps others dealing with these undocumented but useful parts of Matlab.

Allen

[MandarinLu]

A tip: use global function at the beginning for defining the global parameter, so that you can call it in many function handles by just writing global again. For example:

global a
a=10;
%
% here is main code
%
    
function b=Func(location,state)
global a
b=a*100;
end