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

Question

Allen 2020-1-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/498838-use-of-region-and-state-in-pde-models-with-variable-coefficients

评论： MandarinLu 2020-10-27

"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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Ravi Kumar 2020-1-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/498838-use-of-region-and-state-in-pde-models-with-variable-coefficients#answer_408595

在 MATLAB Online 中打开

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

7 个评论
显示 5更早的评论隐藏 5更早的评论

Allen 2020-1-4

在 MATLAB Online 中打开

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?

Allen 2020-1-5

在 MATLAB Online 中打开

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 2020-1-5

在 MATLAB Online 中打开

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 2020-1-6

Ravi,

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

Last line of code should be kOut(state.u<crackT) = Nu*kOut(state.u<crackT);
For some reason if Nu is not defined in kFunc, then the "2 input arguments" error happens.
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.
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 2020-10-27

在 MATLAB Online 中打开

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

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

7 个评论
显示 5更早的评论隐藏 5更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

7 个评论 显示 5更早的评论隐藏 5更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

WeChat

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

7 个评论
显示 5更早的评论隐藏 5更早的评论