Difference between revisions of "Heat equation"

From Medusa: Coordinate Free Mehless Method implementation
Jump to: navigation, search
 
(21 intermediate revisions by 4 users not shown)
Line 1: Line 1:
 
Go back to [[Medusa#Examples|Examples]].
 
Go back to [[Medusa#Examples|Examples]].
 +
 +
This example shows how explicit method of solving partial differential equations is used within the Medusa library. Firstly, this example compares explicit solutions to analytical, then it shows how it can be used to solve heat transfer equation.
  
 
__FORCETOC__
 
__FORCETOC__
 
==Purpose of Heat equation example==
 
 
This example shows how explicit method of solving partial differential equations is used within the medusa library. Firstly, this example compares explicit solutions to analytical, then it shows how it can be used to solve heat transfer equation.
 
  
 
==Domain==
 
==Domain==
  
In this example, a square domain with a cut-out hole in the middle is used. The hole has radius equal to 1/3 of length of a square. Domain nodes are placed randomly with about dx length distance between them.
+
In this example, a square domain with a cut-out hole in the middle is used. The hole has radius equal to 1/3 of length of a square. Domain nodes are placed randomly with about dx length distance between them. Domain boundary is separated in to 5 areas: up, down, left, right and center. The idea is to have Dirichlet boundary condition on the square (up, down, left and right) and Neumann boundary conditions in the middle (center).  
  
<syntaxhighlight lang="c++" line>
+
==Fuction==
BoxShape<Vec2d> box(-3.0, 3.0);
 
DomainDiscretization<Vec2d> domain=box.discretizeBoundaryWithStep(dx);
 
BallShape<Vec2d> ball(0,1.0);
 
DomainDiscretization<Vec2d> domain_empty=ball.discretizeBoundaryWithStep(dx,-5);
 
domain -=domain_empty;
 
GeneralFill<Vec2d> fill_engine;
 
fill_engine(domain, dx);
 
</syntaxhighlight>
 
 
 
Domain boundary is separated in to 5 areas: up, down, left, right and center.
 
 
 
<syntaxhighlight lang="c++" line>
 
auto internal = domain.interior();
 
auto up = domain.types().filter([](int i) { return i == -4 ; });
 
auto down = domain.types().filter([](int i) { return i == -3 ; });
 
auto right = domain.types().filter([](int i) { return  i == -2; });
 
auto left = domain.types().filter([](int i) { return  i == -1; });
 
auto center=domain.types().filter([](int i){ return i == -5;});
 
</syntaxhighlight>
 
 
 
The idea of this domain is to have Dirichlet boundary condition on the square. Furthermore, they are constant over the course of simulation. In the center, the boundary conditions are Neumann.
 
 
 
==Fuction, solution and comparison==
 
  
 
The function used for this example is $u(x,y) = \sin(\pi x)\sin(\pi y)$. Laplace operator of this function is $\Delta u(x,y)=-2\pi^2\sin(\pi x)\sin(\pi y)$.  
 
The function used for this example is $u(x,y) = \sin(\pi x)\sin(\pi y)$. Laplace operator of this function is $\Delta u(x,y)=-2\pi^2\sin(\pi x)\sin(\pi y)$.  
For comparing explicit solution to analytical, we solve function for each node and save solution to u1_analytic. With this we also calculate Dirichlet boundary conditions. Then we set interior and center of domain to arbitrary value, in this case 0. Then we proceed to explicitly solve the domain interior and center. In the domain interior we are solving the following equation:
+
For comparing explicit solution to analytical, we solve function for each node and save solution to $u\_analytic$. With this we also calculate Dirichlet boundary conditions. Then we set domain interior and nodes with Neumann bondary conditions to arbitrary value, in this case 0. We proceed to explicitly solve the domain. In the domain interior we are solving the following equation:
  
 
<math>
 
<math>
Line 45: Line 20:
 
</math>
 
</math>
  
For domain center we set value of $u(x,y,t+1)$ in a way that it satisfies the following equation:
+
For Neumann bondary conditions we set $u(x,y,t+1)$ in a way that it satisfies the following equation:
  
 
<math>
 
<math>
Line 51: Line 26:
 
</math>
 
</math>
  
While doing this, we also compare previous values u1 to new ones u2. If there are differences bigger than dt*1E-5, we continue the simulation. Otherwise, we stop it and compare values in nodes to analytical ones.
+
After the simulation is finished, we can easily calculate relative error. Values of nodes are saved in a vector $u1$, so we can use norm of $|u1 - u\_analytic|$ to calculate error, then divide it by norm of $u1$ to calculate relative error.
 +
 
 +
 
 +
<math>
 +
\begin{align*}
 +
          relative\_error = \frac{\sqrt{\sum{|u1-u\_analytic|^2}}}{\sqrt{\sum{|u1|^2}}}
 +
          \end{align*}
 +
</math>
 +
 
 +
==Weighted Least Squares (WLS)==
 +
 
 +
In this simulation, we use Gaussians for our WLS. Because of scattred nodes in our domain, we have to use 9 basis. The support size we chose is 27. For weight we have chosen 1.0 and for Gaussian weight (sigma) 75. Later we will look how change of sigma and weight effects relative error.
  
 
<syntaxhighlight lang="c++" line>
 
<syntaxhighlight lang="c++" line>
while (ending){
+
WLS<Gaussians<Vec2d>, GaussianWeight<Vec2d>, ScaleToClosest> wls({9, 75.0}, 1.0);
        if (iteracije%izpis==0) {
 
            std::cout << iteracije << std::endl;
 
        }
 
        iteracije++;
 
#pragma omp parallel for
 
        for (int j = 0; j < internal.size(); ++j) {
 
            int i = internal[j];
 
            u2[i] = dt *( op.lap(u1, i) +2*PI*PI*function(domain.pos(i,0),domain.pos(i, 1)))+ u1[i];
 
        }
 
#pragma omp parallel for
 
        for (int j = 0; j < center.size(); ++j) {
 
            int i = center[j];
 
            auto vec=domain.normal(i);
 
            double x=domain.pos(i,0);
 
            double y=domain.pos(i, 1);
 
            u2[i]=op.neumann(u1,i,vec,vec[0]*PI*cos(PI*x)*sin(PI*y)+vec[1]*PI*cos(PI*y)*sin(PI*x));
 
        }
 
        for (int i:internal+center){
 
            double razlika=abs(u1[i]-u2[i]);
 
            if (razlika>konec){
 
                break;
 
            }
 
            ending = false;
 
        }
 
        u1.swap(u2);
 
    }
 
 
</syntaxhighlight>
 
</syntaxhighlight>
  
We ran simulation a couple of times using different dx values. Changing dx value changed number of total nodes as well. We made logarithmic graph of relative error per number of nodes.
+
==Effect of sigma and weight==
 +
 
 +
 
 +
[[File:different_sigma_1.png|600px]]
 +
[[File:different_sigma_3.png|600px]]
 +
 
 +
 
 +
[[File:different_weights_1.png|600px]]
 +
[[File:different_weights_2.png|600px]]
 +
 
 +
Now we use sigma 75 and weight 1.0 to solve heat transfer problem. In this case, we don't have analytical sollution. For our boundary condition, we will set up to 0.0, down to 1.0, right to 0.35 and left to 0.75. For the center hole we will use Neumann boundary conditions. We are solving following equations:
 +
 
 +
<math>
 +
\begin{align*}
 +
        \frac{\partial u (x,y,t)}{\partial t} = \Delta u(x,y,t)
 +
\end{align*}
 +
</math>
 +
 
 +
For the domain interior and for domain center we use:
 +
 
 +
<math>
 +
\begin{align*}
 +
        \frac{\partial u (x,y,t)}{\partial n} = 0
 +
\end{align*}
 +
</math>
  
[[File:error_total_nodes.png|800px]]
+
We get the following picture.
  
Go back to [[Medusa#Examples|Examples]].
+
[[File:heat_trasfer.png|600px]]
 +
[[File:heat_trasfer_2.png|600px]]

Latest revision as of 19:12, 22 October 2022

Go back to Examples.

This example shows how explicit method of solving partial differential equations is used within the Medusa library. Firstly, this example compares explicit solutions to analytical, then it shows how it can be used to solve heat transfer equation.


Domain

In this example, a square domain with a cut-out hole in the middle is used. The hole has radius equal to 1/3 of length of a square. Domain nodes are placed randomly with about dx length distance between them. Domain boundary is separated in to 5 areas: up, down, left, right and center. The idea is to have Dirichlet boundary condition on the square (up, down, left and right) and Neumann boundary conditions in the middle (center).

Fuction

The function used for this example is $u(x,y) = \sin(\pi x)\sin(\pi y)$. Laplace operator of this function is $\Delta u(x,y)=-2\pi^2\sin(\pi x)\sin(\pi y)$. For comparing explicit solution to analytical, we solve function for each node and save solution to $u\_analytic$. With this we also calculate Dirichlet boundary conditions. Then we set domain interior and nodes with Neumann bondary conditions to arbitrary value, in this case 0. We proceed to explicitly solve the domain. In the domain interior we are solving the following equation\[ \begin{align*} u (x,y,t+1) = dt*\bigg(\Delta u(x,y,t) + 2\pi^2\sin(\pi x)\sin(\pi y)\bigg)+u(x,y,t) \end{align*} \]

For Neumann bondary conditions we set $u(x,y,t+1)$ in a way that it satisfies the following equation\[ \frac{\partial u(x,y,t+1)}{\partial n}=n_x \pi \cos(\pi x) \sin(\pi y)+n_y \pi \sin(\pi x) \cos(\pi y) \]

After the simulation is finished, we can easily calculate relative error. Values of nodes are saved in a vector $u1$, so we can use norm of $|u1 - u\_analytic|$ to calculate error, then divide it by norm of $u1$ to calculate relative error.


\( \begin{align*} relative\_error = \frac{\sqrt{\sum{|u1-u\_analytic|^2}}}{\sqrt{\sum{|u1|^2}}} \end{align*} \)

Weighted Least Squares (WLS)

In this simulation, we use Gaussians for our WLS. Because of scattred nodes in our domain, we have to use 9 basis. The support size we chose is 27. For weight we have chosen 1.0 and for Gaussian weight (sigma) 75. Later we will look how change of sigma and weight effects relative error.

1 WLS<Gaussians<Vec2d>, GaussianWeight<Vec2d>, ScaleToClosest> wls({9, 75.0}, 1.0);

Effect of sigma and weight

Different sigma 1.png Different sigma 3.png


Different weights 1.png Different weights 2.png

Now we use sigma 75 and weight 1.0 to solve heat transfer problem. In this case, we don't have analytical sollution. For our boundary condition, we will set up to 0.0, down to 1.0, right to 0.35 and left to 0.75. For the center hole we will use Neumann boundary conditions. We are solving following equations\[ \begin{align*} \frac{\partial u (x,y,t)}{\partial t} = \Delta u(x,y,t) \end{align*} \]

For the domain interior and for domain center we use\[ \begin{align*} \frac{\partial u (x,y,t)}{\partial n} = 0 \end{align*} \]

We get the following picture.

Heat trasfer.png Heat trasfer 2.png