Set a noncontinous 2d heat source in PDE toolbox

hello,

i am trying to set a 2d heatsource in my simulation.

a bit more detail about the problem: i will solve a transient conduction problem of a rotation symmetric cylinder. The source terms are induced by radiation, and come from an external file, created by a different script. The source comes from a 2d matrix, that i want to map onto the PDE toolbox. however i cannot declare a mapping function for qFunc, as i get this error:

Error using pde.HeatSourceAssignment/checkCoefFcnHdlArgCounts (line 240)
Function specifying a heat source must accept two
input arguments and return one output argument.

Error in pde.HeatSourceAssignment/checkFcnHdlArgCounts (line 135)
            self.checkCoefFcnHdlArgCounts(self.HeatSource, systemsize, ndims);

Error in pde.HeatSourceAssignment (line 62)
        self.checkFcnHdlArgCounts(systemsize, numdims);

Error in pde.ThermalModel/internalHeatSource (line 99)
    hs = pde.HeatSourceAssignment(hsContainer,argsToPass{:});

Error in Glass_melt_conduction_only (line 65)
internalHeatSource(thermalmodel,qFunc);

The problem seems to be, that the function is evaluated in the script, but the only thing that is passed on is @heatgeneration. I did it as follows:

...
kFunc = @(region,state) kFunc*region.x;
cFunc = @(region,state) SpecificHeat*region.x;
qFunc = @(region,state) heatgeneration(region)*region.x;
...
internalHeatSource(thermalmodel,qFunc);
...
%%funtion evaluation
function r = getglobal
global q
global Radius
global Height
global n_ring
global n_height

r.a =q;
r.b= Radius;
r.c= Height;
r.d= n_ring;
r.e= n_height;
end

function [result] = heatgeneration(region)
r=getglobal
x=region.x/r.b*r.d;
y=region.y/r.c*r.e;
q=r.a;
if mod(x,1)==0  %if region.x matches the segmentation
    if mod(y,1)==0 %if region.y matches the segmentation
        heat=q(y,x);
    else %if region.y does not match the segmentation
        y1=y-mod(y,1);
        y2=y1+1;
        q1=q(y1,x);
        q2=q(y2,x);
        heat=(q2-q1)/(y2-y1)*(y-y1) + q1;
    end
else %if region.x does not match the segmentation
    if mod(y,1)==0 %if region.y matches the segmentation
        x1=x-mod(x,1);
        x2=x1+1;
        q1=q(y, x1);
        q2=q(y, x2);
        heat=(q2-q1)/(x2-x1)*(x-x1) + q1;
    else %if region.y does not match the segmentation
        y1=y-mod(y,1);
        y2=y1+1;
        x1=x-mod(x,1);
        x2=x1+1;

          q11=q(y1,x1);
          q12=q(y2,x1);
          q21=q(y1, x2);
          q22=q(y2, x2);
          q1=(q12-q11)/(y2-y1)*(y-y1) + q11;
          q2=(q22-q21)/(y2-y1)*(y-y1) + q21;

          heat=(q2-q1)/(x2-x1)*(x-x1) + q1;
      end
  end
  result=heat;
  end

i tried to use global variables, in order to account for the error message, and even tried to make the whole funtion heatgeneration=qfunc, that did not work either.

Answer 1

[Ravi Kumar]

As the message indicates, the function should take 2 input arguments (which your function does, no problem here). However, your function did not return one output argument. The scenario could be that your function definition did not return an output, it errored during execution. There are a couple of places where it might have errored:

1. In the definition: qFunc = @(region,state) heatgeneration(region)*region.x; the variable region.x would be of size 1 x N, where N can be 1 or a large number, meaning it is a rwo vector. The output of heatgeneration(region)*region.x should be of the same size as region.x, please check this first. That means heatgeneration(region) must return a scalar OR you have to use elementwise multiplication heatgeneration(region).*region.x if heatgeneration(region) returns a row vector of same size as that of region.x.

You can debug your heatgeneration(region) to make sure it works when region.x is a row vector. You can test this easily by creating your own region structure as:

region.x = rand(1,10); region.y = rand(1,10); outHeat = heatgeneration(region);

outHeat must be a scalar or a row vector of size 10, otherwise you have an error in the function.

Regards, Ravi

[O Mueller]

Hello Ravi,

Thank you for your answer! I did the test and my function now passes the test you suggested. However: The error message remains. The Problem seems to be, that not the values are being passed into the pde-toolbox, but only the handle @heatgeneration. This function is not defined in the subfunctions. And even if it was:

qFunc = @(region,state) rand(size(region.x,2),1)*region.x;

still gives me an error:

Error using formGlobalKF2D
The number of rows in the f-coefficient matrix, 1500, is not equal the number of
equations in the system, 1.

i am not even sure anymore if i can assign a specific value to each region.

[Ravi Kumar]

Hmm....Lets get this working on a toy problem.

First using the qFunc as you have defined:

thermalmodel = createpde('thermal','transient');
geometryFromEdges(thermalmodel,@squareg);
generateMesh(thermalmodel);
thermalProperties(thermalmodel,'ThermalConductivity',1,'MassDensity',1,'SpecificHeat',1);
thermalBC(thermalmodel,'Edge',1:4,'Temperature',0);
qFunc = @(region,state) rand(size(region.x,2),1)*region.x;
internalHeatSource(thermalmodel,qFunc);
thermalIC(thermalmodel,0);
R = solve(thermalmodel,linspace(0,1,10))

You will get an error.

Next, lets fix the qFunc to make it return correct sized output:

thermalmodel = createpde('thermal','transient');
geometryFromEdges(thermalmodel,@squareg);
generateMesh(thermalmodel);
thermalProperties(thermalmodel,'ThermalConductivity',1,'MassDensity',1,'SpecificHeat',1);
thermalBC(thermalmodel,'Edge',1:4,'Temperature',0);
qFunc = @(region,state) rand(1,size(region.x,1)).*region.x;
internalHeatSource(thermalmodel,qFunc);
thermalIC(thermalmodel,0);
R = solve(thermalmodel,linspace(0,1,10))

Hope these two version of the code would help you resolve your issue.

If you still can't get the error, I would suggest wiring the function qFunc as a normal function:

function qVal = qFunc(region,state)
% code for computing qVal. 
end

You can then insert breakpoints in the function and debug it.

Regards, Ravi

[O Mueller]

Hey Ravi, the working toy-q-funktion is just one random number times region.x It simply assigns one q-value to all regions. Exactly what i don't want to do. I will try to instead of wiring qFunc through heatgeneration, to do it the other way around.

[Ravi Kumar]

I am missing something,. The working toy qFunc assigns size(region.x,1) heat flux values to size(region.x,1) points at which solver wants the heat flux value:

size(region.x,1) will be 1x1 for the initial few calls, where solver is determining the characteristics of the function, if you continue for few more call you will notice region.x will be of size 1xN, N is the order of your number of elements in the model. You may do this experiment by writing the full function as as mentioned at the end of my last post.

The requirement of qFunc is that it must return ONE the value of heat flux at each point point as defined by region.x and region.y vectors, hence heat flux output size should be 1xnumel(region.x).

Better yet, provide a complete reproduction code and sample data, then we can work from there.

Answer 2

[O Mueller]

%%Solve glass melt heat conduction problem
clc;
close all;
clear all;

%%Input

      Radius = 0.3; % [m]
      Height = 0.0635; % [m]
      n_rings = 10; % Number of rings
      n_heights = 20; % Number of heights for which temperature is calculate (At least 2, bottom and top)
      InitialTemperature = 300; % [K] input temperature
      BottomTemperature = 300; % [K] input bottom temperature
      SideTemperature = 300; % [K] input side temperature
      TimeStep = 60; % [s] Timestep for one iteration
      TransientTime = 3600; % [s] Time to calculate transient solution
      SpecificHeat = 1464; % [J/kg/K] Specific heat capacity
      MassDensity = 1202.5; % [kg/m^3] mass density of the glass melt
      kFunc = 0.173; % [W/m/K] 
      MeshSize = 0.003; % [m] Size of the triangular mesh
      HeatFlux = 867.5*ones(n_rings,1); % [W/m^2] Solar irradiance

%%End input section

%%Your input (now array with ones)
load('powerdensity');
q = PD2d;
%we declare now some global variables
setglobal(q,Radius,Height,n_rings,n_heights)

Create a transient thermal model and include the geometry.

thermalmodel = createpde('thermal','transient');

% Define rectangle with different input sources from the radiosity method
for j = 1:n_rings
xVector(j) = j*Radius/n_rings; % Define points on the glass melt surface
yVector(j) = zeros(length(n_rings));
end
g = decsg([2 3+n_rings xVector Radius 0 0 yVector -Height -Height 0]');

figure
pdegplot(g,'EdgeLabels','on','SubdomainLabels','on')
axis equal

% Convert the decsg format into a geometry object. Include the geometry in
% the model.
geometryFromEdges(thermalmodel,g);

% Define Thermal conductivity and specific heat depending on the radius of
% the rod
% data for k taken from pilon 2014 equation 25
kFunc = @(region,state) kFunc*region.x;
cFunc = @(region,state) SpecificHeat*region.x;
%qFunc = @heatgeneration;

% Specify thermal conductivity,  mass density, and specific heat of the
% material.
thermalProperties(thermalmodel,'ThermalConductivity',kFunc,...
                                'MassDensity',MassDensity,...
                                'SpecificHeat',cFunc);

% Define internal heat source (Dein Input qFunc)
internalHeatSource(thermalmodel,qFunc);

% Define boundary conditions, fixed temperatures at the Bottom and side of
% the pot, the BC in the middle r=0 is by default adiabatic!

% thermalBC(thermalmodel,'Edge',1,'Temperature',BottomTemperature);
% thermalBC(thermalmodel,'Edge',2,'Temperature',SideTemperature);

% Create a mesh with elements no larger than MeshSize
msh = generateMesh(thermalmodel,'Hmax',MeshSize);

% Calculate the transient solution
tlist = 0:TimeStep:TransientTime;
thermalIC(thermalmodel,InitialTemperature);
R = solve(thermalmodel,tlist);

for j = 1:n_heights
    tempdisk(j,:) = interpolateTemperature(R,xVector'-Radius/(2*n_rings),yVector'-(j-1)*Height/(n_heights-1),length(tlist):1:length(tlist));
    Y(j,:) = yVector'-(j-1)*Height/(n_heights-1);
    X(j,:) = xVector'-Radius/(2*n_rings);
end

save('Radius','X')
save('Height','Y')
save('tempdisk','tempdisk')

% Function to read heat out of array
function setglobal(q1,Radius1,Height1,n_rings1,n_heights1)
global q
global Rad
global Heig
global n_ri
global n_heig

q =q1;
Rad = Radius1;
Heig = Height1;
n_ri =n_rings1;
n_heig = n_heights1;
end

function r = getglobal
global q
global Rad
global Heig
global n_ri
global n_heig

r.a =q;
r.b= Rad;
r.c= Heig;
r.d= n_ri;
r.e= n_heig;

end

function [result] = qFunc(region, state) %heatgeneration(region, state)
r=getglobal
for p=1:size(region.x,2)
        x=region.x(p)/r.b*r.d 
        y=region.y(p)/r.c*r.e 
        q=r.a;
        y1=y-mod(y,1);
        y2=y1+1;
        x1=x-mod(x,1);
        x2=x1+1;

              if y1==0
                  y1=1;
              end

              if x1==0
                  x1=1;
              end

          if mod(x,1)==0  %if region.x matches the segmentation
              if mod(y,1)==0 %if region.y matches the segmentation
                  heat(p)=q(y,x);
              else %if region.y does not match the segmentation

                  q1=q(y1,x);
                  q2=q(y2,x);
                  heat(p)=(q2-q1)/(y2-y1)*(y-y1) + q1;
              end
          else %if region.x does not match the segmentation
              if mod(y,1)==0 %if region.y matches the segmentation

                  q1=q(y, x1);
                  q2=q(y, x2);
                  heat(p)=(q2-q1)/(x2-x1)*(x-x1) + q1;
              else %if region.y does not match the segmentation

                  q11=q(y1,x1);
                  q12=q(y2,x1);
                  q21=q(y1, x2);
                  q22=q(y2, x2);
                  q1=(q12-q11)/(y2-y1)*(y-y1) + q11;
                  q2=(q22-q21)/(y2-y1)*(y-y1) + q21;

                  heat(p)=(q2-q1)/(x2-x1)*(x-x1) + q1;
              end
          end

end
result=heat*region.x;
end

[Ravi Kumar]

Thanks. You have many coding errors in your qFunc. Before that first need to change the internalHeatSource specification as:

internalHeatSource(thermalmodel,@qFunc);

as qFunc is not a function_handle by itself in this version of your code.

Then insert a breakpoint on the first line of your qFunc and run the code. Step through to resolve all error. I tried to fix a few, but this is not possible without knowing what is your intent. Here are firs few error you will get:

1. heat(p) is assigned without initializing it. This can be fixed by assigning heat = zeros(size(region.x)), before its first use, beginning of the function.

2. Next you will get error because you try to access q(y,x) for x=y=0. Is did the same thing you have done for case x1=y1=0.

3. Then I got the error because q(y,x) is accessed with negative value of y.

Then I stopped. You need to continue debugging this function till the end.

Also at the end you need .* when you multiple heat and region.x as:

result=heat.*region.x;

You may need to consider a better approach to access your tabulated data. Please refer to the interpolation functions.

Regards, Ravi