Problem with coefficient c in function form in solving 2D-PDE (PDE toolbox)

I'm trying to solve a 2D-PDE with a c coefficent in function form, but i'm having some trouble in using location struct.

here is my function for c coefficient:

______________________________________________

function cmatrix = ccoeffunction(location,~)

n1 = 4;

nr = numel(location.x);

Ds=1.4e-10;

jv=2e-5;

cmatrix = zeros(n1,nr);

cmatrix(1,:) = Ds*ones(1,nr);

cmatrix(2,:) = jv*(location.y)+jv*((location.y).^3)/3;

cmatrix(3,:) = -jv*(location.x)+Ds*((location.x)./(location.y));

cmatrix(4,:) = Ds*ones(1,nr);

end

_______________________________________________

However the location struct passed to the function contains in fields x and y just one zero, so my function doesen't work correctly.

______________________________________________

location =

struct with fields:

x: 0

y: 0

z: 0

subdomain: 1

_______________________________________________

How can i pass to my function the correct location struct???

Thanks

Answer 1

[Ravi Kumar]

Hi Chiara,

Solver will call the function several times to check if it is returing the correct format of coeff and so on, for these probing calls a single point is used. If you continue and get past this initial calls, you will be able to see that solver eventually calls with location of several Gauss points where element matrices are integrated.

Regards,

Ravi

[Chiara Moreschini]

The problem is that using the fuction above in my PDE, the solutions does not change even if i change the consant jv in the function. The c coefficient returned is always the same evaluated in location.x=0 and location.y=0:

__________________________

cmatrix =

1.0e-09 *

0.1400

0

NaN

0.1400

__________________________

[Ravi Kumar]

Can you post complete reproduction code, if you haven't been able to setup working.

Regards,

Ravi

[Chiara Moreschini]

Of course, thank you

This is the main program:

__________________________________________

cs0_b=0.2;

Bs=1e-3;

Ds=1.4e-10;

model_sangue = createpde();

rect=[3;4;0;300;300;0;0;0;100;100];

gd=[rect];

ns = char('rect1');

sf='rect1';

dl=decsg(gd);

pdegplot(dl,'EdgeLabels','on','FaceLabels','on')

xlim([0,300])

axis equal;

geometryFromEdges(model_sangue,dl);

pdegplot(model_sangue,'EdgeLabels','on')

applyBoundaryCondition(model_sangue,'dirichlet','Edge',4,'u',cs0_b);

specifyCoefficients(model_sangue,'m',0,'d',1,'c',@ccoeffunction,'a',0,'f',Bs);

setInitialConditions(model_sangue,cs0_b);

generateMesh(model_sangue);

tlist=0:100;

results_sangue = solvepde(model_sangue,tlist);

% PDE 2D plot

cs_b_nodi=results_sangue.NodalSolution;

pdeplot(model_sangue,'XYData',cs_b_nodi,'ZData',cs_b_nodi,'Mesh','on')

xlabel('z')

ylabel('r')

% evaluation of the solution in some random nodes

figure (2);

plot(tlist,cs_b_nodi(450,:), 'b');

xlabel('time');

ylabel('concentrazione in un nodo')

hold on

plot(tlist,cs_b_nodi(2,:), 'y');

hold on

plot(tlist,cs_b_nodi(900,:), 'm');

legend('nodo 450','nodo 2','nodo 900')

__________________________________________

This is the function for c coefficient:

__________________________________________

function cmatrix = ccoeffunction(location,~)

n1 = 4;

nr = numel(location.x);

Ds=1.4e-10;

jv=2e-5;

cmatrix = zeros(n1,nr);

cmatrix(1,:) = Ds*ones(1,nr);

cmatrix(2,:) = jv*(location.y)+jv*((location.y).^3)/3;

cmatrix(3,:) = -jv*(location.x)+Ds*((location.x)./(location.y));

cmatrix(4,:) = Ds*ones(1,nr);

end

__________________________________________

As you can see from the plot of the solution in some random nodes vs time the solution is the same everywere and does not change even if I change the costant jv in the function.

Thank you for your time

Chiara

[Chiara Moreschini]

If you try to use a costant c coefficient like c=[1000;0;-30;1000] instead of ccoeffunction, the solutions change according to the value used in the vector c. I guess this is a confirmation that there's something wrong with the function and the location struct

[Ravi Kumar]

Hi Chiara,

Few items to consider:

1. Nothing is time varying in your setup. In this case you are better of just solve the steay-state problem by setting 'd' to zero and not using tlist in solvepde.
2. If you insert a breakpoint in the ccoeffunction and set a condition, right click and select first option, to stop onlly when numel(location.x)>1, then you will see you get more then one points in location and your c-matrix is compted accordingly.
3. You seem to have been dealing with very small numbers due the value of Ds, Bs, and the size of your domain. Are your coefficients and geometry in consistent units?

Regards,

Ravi

[Chiara Moreschini]

My d coefficient is nonzero because the physical equation i've been trying to resolve has the term du/dt. If i don't use tlist, then i can't see the plot of solution vs time in the random nodes i've choosen.

The value such as Ds are very small because they are physical coefficient (in this case diffusion coefficient) and i cannot modify them; i've also tryed to use grater values to see what happened, but nothing seemed to change. The geometry is in nm, so it is also a small value, according to the other values.

So you are suggesting that my ccoeffunction work properly and the problem is just in the small values i've used?

Thanks

Chiara

[Ravi Kumar]

Well, this could be it! You say your model in nm, I am thining you geometry then shoulod be 300 nm x100 nm. In which case you need to update your geometry creation to account the correct units:

rect=[3;4;0;[300;300;0;0;0;100;100]*1E-9];

If you do so, and plot the last time step solution, then this is what I get:

pdeplot(model_sangue,'XYData',cs_b_nodi(:,end),'ZData',cs_b_nodi(:,end),'Mesh','on')

Based on time-variation, I also think your tlist needs to scaled appropriately.

Regards,

Ravi

[Chiara Moreschini]

Thank you very much!!