Plot along line from. PDE solution

The code below solves a 1-D, transient heat transfer problem set up as in general PDE format. The solution is plotted in color across the domain from 0 to 0.1 after 10 seconds have elapsed. What is the best way to plot the temperature across the length of this domain at this final time?

Thanks

clear all;

%% Create transient thermal model

thermalmodel = createpde(1);

R1= [3,4,0,0.1,0.1,0,0,0,1,1]';
gd= [R1];

sf= 'R1';
ns = char('R1');
ns = ns';
dl = decsg(gd,sf,ns);

%% Create & plot geometry

geometryFromEdges(thermalmodel,dl);
pdegplot(thermalmodel,"EdgeLabels","on","FaceLabels","on")
xlim([0 0.1])
ylim([-1 1])

% axis equal

%% Generate and plot mesh

generateMesh(thermalmodel)

ans = 
  FEMesh with properties:

             Nodes: [2x361 double]
          Elements: [6x152 double]
    MaxElementSize: 0.0402
    MinElementSize: 0.0201
     MeshGradation: 1.5000
    GeometricOrder: 'quadratic'

figure
pdemesh(thermalmodel)
title("Mesh with Quadratic Triangular Elements")

%% Apply BCs

% Edge 4 is left edge; Edge 2 is right side

applyBoundaryCondition(thermalmodel, "dirichlet",Edge=[4],u=100);
applyBoundaryCondition(thermalmodel, "dirichlet",Edge=[2],u=20);

%% Apply thermal properties [copper]

rho= 8933                   % 

rho = 8933

cp= 385                       % 

cp = 385

rhocp= rho*cp             %

rhocp = 3439205

k= 401                         % W/mK

k = 401

%% Define uniform volumetric heat generation rate

Qgen= 0                      % W/m3                      

Qgen = 0

%% Define coefficients for generic Governing Equation to be solved

m= 0

m = 0

a= 0

a = 0

d= rhocp

d = 3439205

c= [k]

c = 401

f= [Qgen]

f = 0

specifyCoefficients(thermalmodel, m=0, d=rhocp, c=k, a=0, f=f);


%% Apply initial condition

setInitialConditions(thermalmodel, 20);


%% Define time limits

tlist= 0: 1: 10;

thermalresults= solvepde(thermalmodel, tlist);

% Plot results

sol = thermalresults.NodalSolution;

subplot(2,2,1)
pdeplot(thermalmodel,"XYData",sol(:,11), ...
    "Contour","on",... 
    "ColorMap","jet")

[Torsten]

What is the best way to plot the temperature across the length of this domain at this final time?

By using "pdepe" instead of the PDE Toolbox.

If you want to continue simulating in 2d for an 1d problem, this code should help:

https://uk.mathworks.com/matlabcentral/answers/52403-plot-solution-of-pde-toolbox-on-a-line-or-submanifold

I guess by "length" you mean from 0 to 0.1 (thus in x-direction).

[John McGrath]

Hi Torsten

Thanks for your suggestion. I tried the plotAlongLine approach that you pointed me to for my problem. Matlab comes back with a message saying it doesn’t recognize p.

I’m actually interested in 2D problems. I thought that by simplifying to a 1D problem I would learn what I wanted to know to apply to 2D cases. It’s probably best to go directly to the type of 2D problem that I’m interested in to ask for your help for plotting along a defined line. Below I include code for a 2D transient problem. It produces a color-coded temperature field at a specified time. I would like to be able to plot the temperature along a line at the specified time. For example, from (0,3) to (6,3). I would also like for the program to define temperature at a specific location and the specific time. For example, T=T(4,3) within the temperature field produced by my program. Thank you in advance.

john

clear all;

thermalmodel = createpde();

R1= [3,4,0,3,3,0,0,0,10,10]';
R2= [3,4,3,6,6,3,-5,-5,5,5]';
gd= [R1 R2];

sf= 'R1+R2';
ns = char('R1','R2');
ns = ns';
dl = decsg(gd,sf,ns);

geometryFromEdges(thermalmodel,dl);
pdegplot(thermalmodel,"EdgeLabels","on","FaceLabels","on")

xlim([-1 7])
ylim([-6 11])
axis equal

%% Generate and plot mesh

generateMesh(thermalmodel);

figure
pdemesh(thermalmodel)
title("Mesh with Quadratic Triangular Elements")

%% Apply BCs

% Edge 3 is left edge; Edge 5 is right side; edge 8 is middle; Edges 2 & 6 top; Edges 1 & 4 bottom
% Edges 7 & 9 offset vertical edges

applyBoundaryCondition(thermalmodel, "dirichlet",Edge=[3],u=100);
applyBoundaryCondition(thermalmodel, "dirichlet",Edge=[5],u=20);


%% Apply thermal properties in first region

rho1= 1                                   % 

rho1 = 1

cp1= 1                                     % 

cp1 = 1

rhocp1= rho1*cp1                  %

rhocp1 = 1

k1= 1                                       % W/mK

k1 = 1

%% Apply thermal properties in second region

rho2= 10                                 % 

rho2 = 10

cp2= 10                                   % 

cp2 = 10

rhocp2= rho2*cp2                  %

rhocp2 = 100

k2= 10                                     % W/mK

k2 = 10

%% Define uniform volumetric heat generation rate

Qgen= 0                      % W/m3                      [1e12]

Qgen = 0

%% Define coefficients for generic Governing Equation to be solved

f= [Qgen]

f = 0

specifyCoefficients(thermalmodel, m=0, d=rhocp1, c=k1, a=0, f=f, Face=1);

specifyCoefficients(thermalmodel, m=0, d=rhocp2, c=k2, a=0, f=f, Face=2);

coeffs= thermalmodel.EquationCoefficients;
findCoefficients(coeffs, Face=1)

ans = 
  CoefficientAssignment with properties:

    RegionType: 'face'
      RegionID: 1
             m: 0
             d: 1
             c: 1
             a: 0
             f: 0

findCoefficients(coeffs, Face=2)

ans = 
  CoefficientAssignment with properties:

    RegionType: 'face'
      RegionID: 2
             m: 0
             d: 100
             c: 10
             a: 0
             f: 0

%% Apply initial condition

setInitialConditions(thermalmodel, 20);


%% Define time limits

tlist = 0:1:50;

%% Solve 

thermalresults= solvepde(thermalmodel, tlist);

% Plot results

sol = thermalresults.NodalSolution;

subplot(2,2,1)
pdeplot(thermalmodel,"XYData",sol(:,50), ...
    "Contour","on",... 
    "ColorMap","jet")

[Torsten]

Plot solution at time 50 from timelist on line connecting [0,0.5] and [6,0.5] using a resolution of 50 points:

plotAlongLine(thermalresults.Mesh.Nodes, thermalresults.NodalSolution(:,50), [0,.5], [6,.5], 50);

[John McGrath]

Matlab returns this:

Unrecognized function or variable 'plotAlongLine'.

[Torsten]

Did you copy the function "plotAlongLine" from the link I gave you and placed it where MATLAB can reach it ?

Maybe you should try MATLAB Onramp before continuing to get the basics of the new programming language:

https://matlabacademy.mathworks.com/details/matlab-onramp/gettingstarted?ef_id=EAIaIQobChMIm5fVjPfBjAMVipuDBx2dJBFBEAAYASAAEgJjWfD_BwE:G:s&s_kwcid=AL!8664!3!548952974296!p!!g!!mathworks%20tutorials&s_eid=ppc_105772043263&q=mathworks+tutorials&gad_source=1

[John McGrath]

I did what I thought was necessary to copy the function where MatLab can reach it but I must have done something wrong.

I completed the MatLab Onramp class. Thanks for the suggestion. It wa squite useful in general but didn't help me to produce the plot I desire.

[Torsten]

You only have to add the call to "plotAlongLine" and the function itself at the end of your code:

clear all;
thermalmodel = createpde();
R1= [3,4,0,3,3,0,0,0,10,10]';
R2= [3,4,3,6,6,3,-5,-5,5,5]';
gd= [R1 R2];
sf= 'R1+R2';
ns = char('R1','R2');
ns = ns';
dl = decsg(gd,sf,ns);
geometryFromEdges(thermalmodel,dl);
pdegplot(thermalmodel,"EdgeLabels","on","FaceLabels","on")
xlim([-1 7])
ylim([-6 11])
axis equal

%% Generate and plot mesh
generateMesh(thermalmodel);
figure
pdemesh(thermalmodel)
title("Mesh with Quadratic Triangular Elements")

%% Apply BCs
% Edge 3 is left edge; Edge 5 is right side; edge 8 is middle; Edges 2 & 6 top; Edges 1 & 4 bottom
% Edges 7 & 9 offset vertical edges
applyBoundaryCondition(thermalmodel, "dirichlet",Edge=[3],u=100);
applyBoundaryCondition(thermalmodel, "dirichlet",Edge=[5],u=20);
%% Apply thermal properties in first region
rho1= 1                                   % 

rho1 = 1

rho1 = 1

rho1 = 1

cp1= 1                                     % 

cp1 = 1

cp1 = 1

cp1 = 1

rhocp1= rho1*cp1                  %

rhocp1 = 1

rhocp1 = 1

rhocp1 = 1

k1= 1                                       % W/mK

k1 = 1

k1 = 1

k1 = 1

%% Apply thermal properties in second region
rho2= 10                                 % 

rho2 = 10

rho2 = 10

rho2 = 10

cp2= 10                                   % 

cp2 = 10

cp2 = 10

cp2 = 10

rhocp2= rho2*cp2                  %

rhocp2 = 100

rhocp2 = 100

rhocp2 = 100

k2= 10                                     % W/mK

k2 = 10

k2 = 10

k2 = 10

%% Define uniform volumetric heat generation rate
Qgen= 0                      % W/m3                      [1e12]

Qgen = 0

Qgen = 0

Qgen = 0

%% Define coefficients for generic Governing Equation to be solved
f= [Qgen]

f = 0

f = 0

f = 0

specifyCoefficients(thermalmodel, m=0, d=rhocp1, c=k1, a=0, f=f, Face=1);
specifyCoefficients(thermalmodel, m=0, d=rhocp2, c=k2, a=0, f=f, Face=2);
coeffs= thermalmodel.EquationCoefficients;
findCoefficients(coeffs, Face=1)

ans = 
  CoefficientAssignment with properties:

    RegionType: 'face'
      RegionID: 1
             m: 0
             d: 1
             c: 1
             a: 0
             f: 0

findCoefficients(coeffs, Face=2)

ans = 
  CoefficientAssignment with properties:

    RegionType: 'face'
      RegionID: 2
             m: 0
             d: 100
             c: 10
             a: 0
             f: 0

%% Apply initial condition
setInitialConditions(thermalmodel, 20);
%% Define time limits
tlist = 0:1:50;
%% Solve 
thermalresults= solvepde(thermalmodel, tlist);
% Plot results
sol = thermalresults.NodalSolution;
subplot(2,2,1)
pdeplot(thermalmodel,"XYData",sol(:,50), ...
    "Contour","on",... 
    "ColorMap","jet")

plotAlongLine(thermalresults.Mesh.Nodes, thermalresults.NodalSolution(:,50), [0,.5], [6,.5], 50);

function plotAlongLine(p, u, xy1, xy2, numpts)
x = linspace(xy1(1),xy2(1),numpts);
y = linspace(xy1(2),xy2(2),numpts);
F = TriScatteredInterp(p(1,:)', p(2,:)', u);
uxy = F(x,y);
figure; plot(x, uxy);
end

[John McGrath]

I copied what you suggested above and pasted it into the end of my code

plotAlongLine(thermalresults.Mesh.Nodes, thermalresults.NodalSolution(:,50), [0,.5], [6,.5], 50);

function plotAlongLine(p, u, xy1, xy2, numpts)

x = linspace(xy1(1),xy2(1),numpts);

y = linspace(xy1(2),xy2(2),numpts);

F = TriScatteredInterp(p(1,:)', p(2,:)', u);

uxy = F(x,y);

figure; plot(x, uxy);

end

and MATLAB returned this:

Unrecognized function or variable 'plotAlongLine'.

What you suggested makes perfect sense to me but I don't see what I am doing wrong

thanks

[Torsten]

I'm surprised you have these problems. Last year, the function worked fine for you:

https://uk.mathworks.com/matlabcentral/answers/2149574-how-to-plot-temperature-along-a-defined-line-within-a-2d-temperature-field?s_tid=srchtitle

Which MATLAB version do you use ?

If it's an old one, you could try removing the function at the end of your code and instead saving it in your MATLAB working directory as "plotAlongLine.m" .

If this does not help, I don't know. You see that it works here using MATLAB online without problems.

Does this code work on your MATLAB ?

a = 4;
a_squared = fun(a)

asquared = 16

function result = fun(a)
  result = a.^2;
end

[John McGrath]

Last year the way I defined my problem was as an explicit heat transfer problem. That worked until I tried to expand my problem to include time varying convective boundary conditions. When I ran into problems doing that the advice I got from the MATLAB community was that I needed to set up my problem as a general pde problem rather than a heat transfer problem per se.

I write my code in a Word document then copy and paste it into MATLAB.

The code you suggested will not work with my MATLAB [MATLAB_R2024a] if pasted directly into the Command Window. It will work if copied and pasted into a New Script and then Run Section is executed. Same is true if it is copied and pasted into a New Live Script

I have discovered that some code works pasted into the Command Window directly. The same code copied and pasted into the Command WIndow directly will not work in some cases but will work if pasted into the New Script or New Live Script areas. I'm trying to understand why.

[Torsten]

You cannot execute MATLAB code containing self-written functions by pasting it into the command window. The function is not known at the moment it is called.

In Octave, it works if you save a function under its name in the working directory. In this case, it can subsequently be called from the command window.

[John McGrath]

What is Octave?

thanks for your continued support

[Torsten]

What is Octave?

https://octave.org/

Answer 1

[Karanjot]

Hi John,

In response to your latest comment, MATLAB returns the following error,as 'plotAlongLine' is not a standard function that is part of MATLAB. It's a custom function defined by another user.

Unrecognized function or variable 'plotAlongLine'.

In order to resolve the error, you may copy the function definition into your codebase, from the link attached below:

function plotAlongLine(p, u, xy1, xy2, numpts)
x = linspace(xy1(1),xy2(1),numpts);
y = linspace(xy1(2),xy2(2),numpts);
F = TriScatteredInterp(p(1,:)', p(2,:)', u);
uxy = F(x,y);
figure; plot(x, uxy);
end

https://www.mathworks.com/matlabcentral/answers/52403-plot-solution-of-pde-toolbox-on-a-line-or-submanifold#answer_64119

As also suggested by Torsten, Feel free to go through the following beginner courses on MATLAB basics:

https://matlabacademy.mathworks.com/?page=1&fq=get-started-with-matlab&sort=featured

I hope this helps!

[John McGrath]

Hi Karanjot,

Thank you very much for your suggestions. I have completed some of these courses and will continue to work on them.

I tried to implement your suggestion for code- but unsuccessfully. Furthermore, when I searched for documentation regarding the TriScatteredInterp function, the MatLab documentation indicates that this function is not recommended and that the function scatteredInterpolant should be used instead. [The original code taht incldued this functoin dates abck to 2012 or something like that.

I initially defined my heat transfer model in terms of a thermal model. I ran into problems when I tried to implement time-varying convective boundary conditions. I was told that to avoid that issue, I should define my heat transfer problem as a general pde problem- which is what I am trying to do now.

I recognize that I do not understand the details of how my pde solution for the temperatures [1D or 2D] [steady or transient] are mapped onto the grid that is generated by the generateMesh function. I therefore don't have a clear picture of how to query the vectors or matrices that are being generated so as to produce values at defined locations or to produce desired plots.

Thanks again

john