John_3265
John_3265

Reputation: 71

1D Heat Diffusion PDE implementation in Modelica(Dymola)

I was trying to implement a 1D heat diffusion example from the Peter Fritzon's "Object Oriented Modeling and Simulation with Modelica 3.3" based on the Grid Function Finite Difference Approach to 1D Heat Diffusion but I am getting an error:

Translation of Heat_diffusion_Test_1D.HeatDiffusion1D:

The DAE has 50 scalar unknowns and 50 scalar equations.

Cannot find differentiation function:
Heat_diffusion_Test_1D.DifferentialOperators1D.pder(
u, 
1)
with respect to time

    Failed to differentiate the equation
    Heat_diffusion_Test_1D.FieldDomainOperators1D.right(Heat_diffusion_Test_1D.DifferentialOperators1D.pder (
    u, 
    1)) = 0;

    in order to reduce the DAE index.

    Failed to reduce the DAE index.

    Translation aborted.

    ERROR: 1 error was found

    WARNING: 1 warning was issued

I am a quite new user of the programming language, do you have any idea on how I could fix the issue? Apparently the problem is in the boundary condition of the right boundary. Thank you in advance!

here is the code:

package Heat_diffusion_Test_1D
import SI = Modelica.SIunits;
  model HeatDiffusion1D
    import SI = Modelica.SIunits;
    parameter SI.Length L=1;
    import Heat_diffusion_Test_1D.Fields.*;
    import Heat_diffusion_Test_1D.Domains.*;
    import Heat_diffusion_Test_1D.FieldDomainOperators1D.*;
    import Heat_diffusion_Test_1D.DifferentialOperators1D.*;
    constant Real PI=Modelica.Constants.pi;

    Domain1D rod(left=0, right=L, n=50);
    Field1D u(domain=rod, redeclare type FieldValueType=SI.Temperature, start={20*sin(PI/2*x)+ 300 for x in rod.x});

  equation 
    interior(der(u.val))=interior(4*pder(u,x=2));
    left(u.val)=20*sin(PI/12*time)+300;
   right(pder(u,x=1)) = 0;
    //right(u.val)=320;

    annotation (Icon(coordinateSystem(preserveAspectRatio=false)), Diagram(
          coordinateSystem(preserveAspectRatio=false)));
  end HeatDiffusion1D;

  package Fields

    record Field1D "Field over a 1D spatial domain"
      replaceable type FieldValueType = Real;
      parameter Domains.Domain1D domain;
      parameter FieldValueType start[domain.n]=zeros(domain.n);
      FieldValueType val[domain.n](start=start);

      annotation (Icon(coordinateSystem(preserveAspectRatio=false)), Diagram(
            coordinateSystem(preserveAspectRatio=false)));
    end Field1D;
    annotation ();
  end Fields;

  package Domains
    import SI = Modelica.SIunits;
    record Domain1D "1D spatial domain"
      parameter SI.Length left=0.0;
      parameter SI.Length right=1.0;
      parameter Integer n=100;
      parameter SI.Length dx = (right-left)/(n-1);
      parameter SI.Length x[n]=linspace(right,left,n);

      annotation (Icon(coordinateSystem(preserveAspectRatio=false)), Diagram(
            coordinateSystem(preserveAspectRatio=false)));
    end Domain1D;
    annotation ();
  end Domains;

  package FieldDomainOperators1D 
    "Selection operators of field on 1D domain"

    function left "Returns the left value of the field vector v1"
      input Real[:] v1;
      output Real v2;

    algorithm 
      v2 := v1[1];
    end left;

    function right "Returns the left value of the field vector v1"
      input Real[:] v1;
      output Real v2;

    algorithm 
      v2 := v1[end];

    end right;

    function interior 
      "returns the interior values of the field as a vector"
      input Real v1[:];
      output Real v2[size(v1,1)-2]; 

    algorithm 
      v2:= v1[2:end-1];

    end interior;


  end FieldDomainOperators1D;

  package DifferentialOperators1D 
    "Finite difference differential operators"

    function pder 
      "returns vector of spatial derivative values of a 1D field"
      input Heat_diffusion_Test_1D.Fields.Field1D f;
      input Integer x=1 "Diff order - first or second order derivative";
      output Real df[size(f.val,1)];

    algorithm 
      df:= if x == 1 then SecondOrder.diff1(f.val, f.domain.dx)
     else 
         SecondOrder.diff2(f.val, f.domain.dx);

    end pder;

    package SecondOrder 
      "Second order polynomial derivative approximations"

      function diff1 "First derivative"
        input Real u[:];
        input Real dx;
        output Real du[size(u,1)];

      algorithm 
        du := cat( 1, {(-3*u[1] + 4*u[2] - u[3])/2/dx},  (u[3:end] - u[1:end-2])/2/dx, {(3*u[end] -4*u[end-1] + u[end-2])/2/dx});

      end diff1;

      function diff2
        input Real u[:];
        input Real dx;
        output Real du2[size(u,1)];
      algorithm 
        du2 :=cat( 1, {(2*u[1] - 5*u[2] + 4*u[3]- u[4])/dx/dx}, (u[3:end] - 2*u[2:end-1] + u[1:end-2])/dx/dx, {(2*u[end] -5*u[end-1] + 4*u[end-2] - u[end-3])/dx/dx});
      end diff2;
      annotation ();
    end SecondOrder;

    package FirstOrder 
      "First order polynomial derivative approximations"

      function diff1 "First derivative"
        input Real u[:];
        input Real dx;
        output Real du[size(u,1)];

      algorithm 
                // Left, central and right differences
        du := cat(1, {(u[2] - u[1])/dx}, (u[3:end]-u[1:end-2])/2/dx, {(u[end] - u[end-1])/dx});

      end diff1;
      annotation ();
    end FirstOrder;
    annotation ();
  end DifferentialOperators1D;
  annotation (uses(Modelica(version="3.2.1")));
end Heat_diffusion_Test_1D;

Upvotes: 2

Views: 533

Answers (1)

Lukas Exel
Lukas Exel

Reputation: 86

Your implementation translate and simulate in Dymola 2015 FD01 if you inline all functions in the package FieldDomainOperators1D and DifferentialOperators1D with annotation(Inline=true). The reason is given by Peter Fritzson at page 406:

[...] Note that inline semantics is needed for the pder() function and the grid function. [...] This makes it possible to put calls to such a function at places where function calls are usually not allowed. [...]

Upvotes: 3

Related Questions