Does hyperbolic-solver really use a nonlinear solver?!

My PDE is time-dependent and has a coefficient c , which is depending on the gradient of the solution u (due to the usage of Green-Lagrange-strain measure).

Following the documentation „a nonlinear solver is available for the nonlinear elliptic PDE, the parabolic and hyperbolic equation solvers also solve nonlinear and time-dependent problems.“ For elliptic problems it is necessary to explicitly use pdenonlin instead of assempde . But according to the documentation „the parabolic and hyperbolic functions call the nonlinear solver automatically “.

I don't use the PDE-Tool-GUI, but implemented my own method of line with the pdetool-command-line functions. I looked up the codes of the hyperbolic / parabolic -solver and it's subroutines, and I can nowhere find any call of a nonlinear-solver in case of solution-dependent coefficients!

The results of my calculations also strongly imply, that the nonlinearity induced by c is not treated correctly, since the results are the same as in the linear case. So if anyone can tell me, if it is possible to solve nonlinear transient PDEs with the pdetool-command-line functions, I would be very happy, since I am at my wit's end for another error source producing the unexpected results.

Answer 1

[Bill Greene]

Hi,

Unfortunately, I believe you have uncovered a bug in the hyperbolic function that occurs when you have a nonlinear equation but no Dirichlet boundary conditions. I believe there is a simple fix. If you replace the line

uFull = u;

in toolbox/pde/pdehypdf.m with

uFull = u(1:nu);

your example will run to completion. The line to replace is line 16 in MATLAB R2013b.

However, the results may not be what you intended. I'm guessing that your problem is just a rigid body rotation of the block. If so, the rotated block is larger than the original block which is obviously not correct for the rigid body case.

Before saying more, I should first say that I have wanted to derive a c-matrix for a Lagrangian formulation of 2D elasticity but have been somewhat discouraged by all the algebra required. So I don't have such formulation but believe it is possible to derive one and suspect it is more complicated than what you have in cCoeff_rotation.

I did a simple test of your cCoeff_rotation function (file attached) where I created a simple block, defined the displacement vector for a 45-degree rotation, and calculated the internal forces using cCoeff_rotation. This vector is not zero as expected.

You probably have already looked at this page:

http://www.mathworks.com/help/pde/ug/nonlinear-equations.html?searchHighlight=nonlinear+equations

but there is a key equation there that is easy to overlook. Under the line "The residual vector ρ(U) can be easily computed as" you will see that the residual vector is calculated from a product of K and U. This is different from the typical nonlinear mechanics notation where K is a tangent stiffness matrix.

Hope this documentation page and the fix will point you in the right direction.

Regards,

Bill

[Lena]

Thank you for this good example, it really filters out the main problem! So I checked my calculation of c one more time, and you were right: There was really an error in the calculation. I've attached the corrected version. The internal force has descreased a little bit, but unfortunately is still far away from zero.

By the way there seems to be a inconsistency in the notation on page 2-51 in the Dokumenation : One the one hand we have c(2N+1) in 3rd column and on the other hand we have c(5) in the 3rd line of the Tensor again.

Do you think there could perhaps be an error in my referring of the gradient? Thanks for your great help, Lena

[Lena]

The complicated Notation on page 2-50 is to understand as a product rule, isn't it? E.g. the first summand is

d/dx*c(ij11)*duj/dx =

= d/dx (c(ij11) duj/dx) = dc(ij11)/dx*duj/dx + c(ij11)*d2uj/dx2

If not so, I have to repeat my calculations for c....

[Lena]

I have just found the error!!! Instead of

u = u(:);

We of course have to write

u = [u(1,:)'; u(2,:)'];

Now we get

Max internal force = 3.05176e-05

I'm now trying to solve my original Problem, and inform you about my results! Thank you Bill!

[Hassan Diab]

I am using matlab 2011, I have the pdenonlin command but i didn't find the line to change,

I have the same problem, pdenonlin try to divide H and get a memory error when there is not any Dirichlet condition.

do you think this issue is solved in the new release of Matlab.

Thanks

Answer 2

[Bill Greene]

The first thing to check is that your version of MATLAB is R2012b or newer. The nonlinear hyperbolic (and parabolic) solvers have been available only since then.

Since you didn't post your code for computation of the c-coefficient, I'm speculating about another possibility. The hyperbolic (and parabolic) functions detect that c is a function of u (or ux, uy) by passing a vector, u=NaN, to your routine for calculating c. They expect that the computed c-matrix contains at least one NaN. Usually this requirement is satisfied automatically because computations involving NaN produce a NaN as a result. But occasionally you need to treat this NaN call specially.

Finally, I suggest you place a debugger break point in your routine that calculates the c-coefficient and then run your example problem. Obviously, the hyperbolic solver should be calling this routine many times during the analysis and each time it will stop at your break point.

If none of these pointers help you resolve the issue, please post your code and I'll take a look.

Bill

[Lena]

Thanks for your quick answer! I've checked the things you recommended, but the difficulties continue to exist. As a model problem I considered a rotation around one fixed point, so we consider large deformations. In the attached example I left this fixed point out for simplicity, since otherwise it has to be given manually to the odesolver. But the main effect is also visible from this non-fixed case.

I don't use the hyperbolic solver itself, because I want to implement my own time integrator. Only for this example I used the ode15s, which I have replaced in my work by another solver. But the problem of giving over the right-hand-side, mass and jacobi is the same. I have just proceeded as in the code of hyperbolic.m and it's subroutines.

If you run the code with both c = @cCoeff_rotation; and c = @cCoeff_rotation_linear; you get the same results: The block moves down and it's size increases!

This increasing of the size under rotation should be avoided by use of the nonlinear coefficient c (nonlinear strain measure), so I imply there was never used a nonlinear solver.

function [] = rotation()
clear all; close all;
global p; global t; global e; global c; global b; global a; global f;
global E1; global nu1; global Q0; global K0;

Tspan = [0.0, 6];
E1 = 7.0e10;
rho1 = 2700;
nu1 = 0.34;

%Coefficients
a = 0.0;
d = rho1;
f = [0.0; 0.0];
c = @cCoeff_rotation;
%c = @cCoeff_rotation_linear;
b = @wbound_rotation;

%Geometry
gd = [2.0000;5.0000;-1.0000;-0.5000;1.0000;1.0000;-1.0000;-0.5000;...
    -0.5000;-0.5000;0.5000;0.5000];
ns =[ 80; 49];
sf = 'P1';
gd = decsg(gd, sf, ns);

%Mesh
p =[1.0000 1.0000 -1.0000 -1.0000 -0.5000;
    -0.5000 0.5000 0.5000 -0.5000 -0.5000];
e = [1 2 3 4 5; 2 3 4 5 1; 0 0 0 0 0; 1 1 1 1 1; 1 2 3 4 5; 1 1 1 1 1;...
    0 0 0 0 0];
t = [1 2 3; 2 3 4; 5 5 5; 1 1 1];

np = size(p,2);
ne = size(e, 2);
u0 = zeros(2*np, 1);
ut0 = zeros(2*np, 1);

%Space-discretization
t0 = min(Tspan);
[~, MM0] = assema(p, t, 0, d, f, t0);
[K0, M0, F0] = assema(p, t, c, a, f, u0, t0);
[Q0, G0, ~, ~] = assemb(b, p, e, u0, t0);
nu = size(MM0,1);
mass = [speye(nu,nu) sparse(nu,nu); sparse(nu,nu) MM0];

%Odesovler
uu0 = [u0; ut0];
options = odeset('Jacobian', @Jac_rotation);
options = odeset(options, 'Mass', mass);
reltol = 1e-3;
options = odeset(options, 'RelTol', reltol, 'Stats', 'on', 'MaxOrder', 2);
[tt, uu] = ode15s(@ODEfct_rotation, Tspan, uu0, options);
uu = uu';
uSol = uu(1:2*np,:);

%Plots
result = uSol(:,end);
pdeplot(p,e,t)
hold on
pdeplot(p+[result(1:np) result(np+1:np+np)]', e, t, 'mesh', 'on', 'title', 'deformed mesh at last timestep');
xlabel('x-axis')
ylabel('y-axis')
axis equal
end

function [c] = cCoeff_rotation(p, t, u, time)
global E1;
global nu1;
[UX, UY] = pdegrad(p,t,u);
u1x = UX(1,:);
u1y = UY(1,:);
u2x = UX(2,:);
u2y = UY(2,:);
g1 = (1-nu1)/2;
e1 = E1/(1-nu1^2);
tn = size(t, 2);
c = zeros(16, tn);
for i=1:tn
    if(isnan(time) || any(isnan(u)))
        c(:, i) = NaN*ones(16, 1);
    else
        c(:, i) = e1*[1+0.5*u1x(i);
            0.5*nu1*u1y(i);
            g1*u1y(i);
            g1;
            0.5*u2x(i);
            nu1+0.5*nu1*u2y(i);
            g1+g1*u2y(i);
            0;
            g1*u1y(i);
            g1;
            nu1+0.5*nu1*u1x(i);
            0.5*u1y(i);
            g1+g1*u2y(i);
            0;
            0.5*nu1*u2x(i);
            1+0.5*u2y(i)];
    end
end
end

function [c] = cCoeff_rotation_linear(p, t, u, time)
global E1;
global nu1;
Gm1 = E1/(2*(1+nu1));
mu1 = 2*Gm1*(nu1/(1-nu1));
tn = size(t, 2);
c = zeros(16, tn);
for i=1:tn
    c(:,i) = [2*Gm1+mu1 0 0 Gm1, 0 mu1 Gm1 0, 0 Gm1 mu1 0, Gm1 0 0 2*Gm1+mu1]';
end
end

function [df] = Jac_rotation(time,u) % according to pdehypdf.m
global Q0; global p;  global c; global a; global t;
nu = size(Q0,1);
uu =  u(1:nu);
[K0, M0, ~] = assema(p, t, c, a, [0;0], uu, time);
KK = K0+M0+Q0; %only G0 varies due to time-dependent NBC
df = - [sparse(nu,nu) -speye(nu,nu); KK sparse(nu,nu)];
end

function [fct]= ODEfct_rotation(time,u) % according to pdehypf.m
global K0; global b; global p;  global e; global c; global a;
global f; global t;
nu = size(K0,1);
uu =  u(1:nu);
[K0, M0, F0] = assema(p, t, c, a, f, uu, time);
[Q0, G0, ~, ~] = assemb(b, p, e, uu, time);
KK = K0 + M0 + Q0;
K = [sparse(nu,nu) -speye(nu,nu); KK sparse(nu,nu)];
FF =F0+G0;
F = [sparse(nu,1); FF];
fct = -K*u+F;
end

function [q,g,h,r]=wbound_rotation(p,e,u,time)
%WBOUND_ROTATION  Boundary condition data.
bl=[
    2 2 2 2 2
    0 0 0 0 0
    1 1 1 1 1
    1 1 1 1 1
    1 1 1 1 1
    1 1 1 1 1
    1 1 1 20 1
    1 1 1 20 1
    48 48 48 48 48
    48 48 48 48 48
    48 48 48 48 48
    48 48 48 48 48
    48 48 48 40 48
    48 48 48 53 48
    48 48 48 48 48
    48 48 48 48 48
    48 48 48 42 48
    48 48 48 115 48
    48 48 48 105 48
    48 48 48 110 48
    49 49 49 40 49
    48 48 48 48 48
    48 48 48 46 48
    49 49 49 53 49
    48 48 48 42 48
    48 48 48 116 48
    0 0 0 41 0
    0 0 0 41 0
    0 0 0 46 0
    0 0 0 42 0
    0 0 0 110 0
    0 0 0 120 0
    0 0 0 40 0
    0 0 0 53 0
    0 0 0 48 0
    0 0 0 48 0
    0 0 0 42 0
    0 0 0 115 0
    0 0 0 105 0
    0 0 0 110 0
    0 0 0 40 0
    0 0 0 48 0
    0 0 0 46 0
    0 0 0 53 0
    0 0 0 42 0
    0 0 0 116 0
    0 0 0 41 0
    0 0 0 41 0
    0 0 0 46 0
    0 0 0 42 0
    0 0 0 110 0
    0 0 0 121 0
    ];
if any(size(u))
    [q,g,h,r]=pdeexpd(p,e,u,time,bl);
else
    [q,g,h,r]=pdeexpd(p,e,time,bl);
end
end

Answer 3

[Bill Greene]

Oh, sorry, I thought you were using the documented hyperbolic function. I'm afraid I have no idea why your code isn't doing what you expect.

Bill

[Lena]

I have just tried to write my problem in the pdetool-GUI (=> use hyperbolic function). So I have to enter the coefficiets c11 , c12 ,... I tried for example to refer to the first component du1/dx of the gradient with ux(1,:) e.g.

c11 = 7.9150e10.*(1+0.5.*ux(1,:)) 7.9150e10.*(0.5*0.34.*uy(1,:)) 7.9150e10.*(0.33.*uy(1,:)) 7.9150e10*0.33

But I got "PDE Toolbox Error: Matrix dimensions must agree". What is the correct handling for the gradient in the system case? I think that, this is also one problem in my code.

Another remark to my code: My implementation follows the code hyperbolic.m . In the code of hyperbolic.m and it's subroutines I can nowhere find any call of a nonlinear solver such as pdenonlin.m , which uses Gauss-Newton or something like that to treat the nonlinear equation. Where exactly in hyperbolic.m (or parabolic.m ) is treated a nonlinearity like K(u)u=F like in the elliptic case? Thank you very much, Lena

Answer 4

[Bill Greene]

Hi,

I think that referring to gradient components in the form e.g. ux(1,:) is fine. But entering complicated nonlinear coefficients for the system (N=2) case in the pdetool GUI can be tricky. I suggest solving the problem using command-line functions instead of the GUI.

I have attached code for a simple linear cantilever beam transient response example that uses the hyperbolic function. If you already have expressions for the geometrically nonlinear c-coeffients, you have already done most of the work. It should be straightforward to code those in the calcCmat function in the example.

Bill

[Lena]

Thank you, where can I get the code?

Answer 5

[Bill Greene]

Sorry about that! Apparently this GUI is too sophisticated for me ;-) I think I got it this time.

Bill

[Lena]

Hi Bill, thank you for your code, I have adapted my code to the hyperbolic-solver now (but keeping in mind that I want to replace its ode-solver). I simply used my implementation of the c -Coefficient from the above code:

Line 15: c = @cCoeff_rotation;
Line 16: %c = @cCoeff_rotation_linear;

If you run it with Line 15 (the nonlinear Version) I get the error message

Error using  + 
Matrix dimensions must agree.

With Line 16 it runs, but with the block's size increasing due to the use of a linear strain measure.

What is wrong in the nonlinear c-Coefficient-Function? Thank you very much for your help! Lena