Skip to content
Snippets Groups Projects
Commit aef0d678 authored by Lars Hubatsch's avatar Lars Hubatsch
Browse files

Matlab and Python agree. Fixed bug in power notation.

parent c308a6ee
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
@Ternary_model/plot_sim.m
+ 1
− 1
View file @ aef0d678
......@@ -24,7 +24,7 @@ if strcmp(mode, 'mov')
elseif strcmp(mode, 'plot')
% figure;
hold on;
xlim([4, 8]); ylim([0, 2]);
xlim([0, 4]); ylim([0, 2]);
ax = gca;
ax.FontSize = 26;
xlabel('x [\mum]'); ylabel('volume fraction');
......
Like1D_Python_Matlab.ipynb
+ 41
− 106
View file @ aef0d678
%% Cell type:code id: tags:
``` python
from __future__ import print_function
from fenics import *
from mshr import *
from dolfin import *
import matplotlib.pyplot as plt
import time
from math import *
import numpy as np
from scipy.interpolate import interp1d
def Load_matlab():
### Import Matlab Results ###
Matlab_result = np.loadtxt('resultsmatlab.txt')
Matlab_mesh = np.loadtxt('meshematlab.txt')
def Solve():
set_log_level(40)
dt = 0.1
##PARAMETERS##
nu = 10E-6
nu = 10**-6
Xi = 7./3
a = -1.
b = sqrt(Xi/(Xi -2))*pow(nu,1/3)
e = sqrt((3./8)*(Xi - 2))
u0 = 0.5
e_g0 = 1.
u_g0 = 0.5
### Import Matlab Results ###
Matlab_result = np.loadtxt('resultsmatlab.txt')
Matlab_mesh = np.loadtxt('meshematlab.txt')
Matlab_result = genfromtxt('resultsmatlab.txt', delimiter=',')
Matlab_mesh = genfromtxt('meshematlab.txt', delimiter=',')
Matlab_phi_t = genfromtxt('phi_tot.txt', delimiter=',')
############ Create The Problem ########
#mesh = Mesh('Mesh1D.xml')
mesh = Mesh('1DMesh_Larsidea4.xml')
mesh = Mesh('Meshes/Like1DMesh.xml')
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh,P1)
Phi_u = Function(V)
v = TestFunction(V)
Phi_u0 = Function(V)
Phi_tot = Function(V)
Phi_tot = interpolate(Expression('e*tanh(-(x[0]+a)/b)+u0', degree=2,e = e, a = a, b = b, u0 = u0),V)
Phi_u0 = interpolate(Expression('x[0] <= -a ? 0.0 : (u0 - e)' , degree=2, a = a, e = e, u0 = u0),V)
Gamma_0 = interpolate(Expression('e_g0*(tanh((x[0]+a)/b)+1)/2+u_g0;',degree = 2, e_g0 = e_g0, a = a, b = b,u_g0 = u_g0), V)
####### CREATE MATRIX FOR PLOT THE RESULTS FROM PYTHON AT EACH TIME ################
tol = .001 # avoid hitting points outside the domain
x = np.linspace(0 + tol, 4 - tol, 20000)
points = [(x_, 0.01) for x_ in x]
Matrix_result_Python = np.array([Phi_u0(point) for point in points])
#### PLOT IC MATLAB
Value_Matlab = []
X_Matlab = []
for j in range(1,4124):
Value_Matlab.append(Matlab_result[0][j])
X_Matlab.append(Matlab_mesh[j])
plt.plot(X_Matlab,Value_Matlab, label = "Matlab")
#### PLOT IC PYTHON
Value_Python =[]
X_Python = []
tol = .001 # avoid hitting points outside the domain
x = np.linspace(0.99 + tol, 1.01 - tol, 500000)
points = [(x_, 0.01) for x_ in x] # 2D points
Phi_u0_x1 = np.array([Phi_u0(point) for point in points])
plt.plot(x, Phi_u0_x1, label = "Python")
### PLOT IC ###
plt.axis([0.9,1.1,-0.05,0.2])
plt.ylabel('Phi_u0')
plt.xlabel('x')
# Gamma_0 = interpolate(Expression('e_g0*(tanh((x[0]+a)/b)+1)/2+u_g0;',degree = 2, e_g0 = e_g0, a = a, b = b,u_g0 = u_g0), V)
#### PLOT IC MATLAB AND PYTHON
x = np.linspace(0, 4, 10000)
plt.plot(x, [Phi_tot(x, 0.01) for x in x], label = 'Python')
plt.plot(Matlab_mesh, Matlab_phi_t, label = 'Matlab')
plt.ylim(0, 1)
plt.xlim(0, 4)
plt.legend()
plt.show()
### Boundary Conditions ###
u_D = Expression('u0 - e', degree = 2, u0 = u0, e=e)
def boundary_right(x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0],30,tol)
bc = DirichletBC(V, u_D, boundary_right)
### Form ###
Form = (inner((Phi_u-Phi_u0)/dt, v) - Gamma_0*(inner(((1-Phi_tot)/Phi_tot)*Phi_u*grad(Phi_tot) - (1-Phi_tot)*grad(Phi_u), grad(v)))) *dx
Form = (dot((Phi_u-Phi_u0)/dt, v)*dx -
(1-Phi_tot)/Phi_tot*Phi_u*dot(grad(Phi_tot), grad(v))*dx +
(1-Phi_tot)*dot(grad(Phi_u), grad(v))*dx)
### Stock Phiu0, Phitot ###
vtkfile = File('like_1D/Phi0.pvd')
vtkfile << Phi_u0
vtkfile = File('like_1D/Phitot.pvd')
vtkfile << Phi_tot
#### Open the file XMDF for writing results ####
t = 0
cFile = XDMFFile('like_1D.xmdf')
cFile.write(Phi_u0, t)
ti = time.time()
### Computing ###
for i in range(10):
for i in range(100):
### plot results ###
if i % 10 == 0:
tol = .001 # avoid hitting points outside the domain
x = np.linspace(0 + tol, 4 - tol, 20000)
points = [(x_, 0.01) for x_ in x]
Phi_u0_x1 = np.array([Phi_u0(point) for point in points])
plt.plot(x, Phi_u0_x1, label = "Python")
plt.plot(Matlab_mesh, Matlab_phi_t, label = 'Phi_tot')
plt.plot(Matlab_mesh, Matlab_result[int(i), :], linestyle='--', lw=2,
label = 'Matlab')
### Plot ###
plt.axis([0,2,0.0,0.4])
plt.ylabel('Phi_u0')
plt.xlabel('x')
plt.legend()
plt.show()
solve(Form == 0, Phi_u, bc,solver_parameters={'newton_solver': {'linear_solver': 'gmres'}})
assign(Phi_u0, Phi_u)
t += dt
cFile.write(Phi_u0, t)
print(time.time() - ti)
### plot Matlab Results ###
Value_Matlab = []
X_Matlab = []
for j in range(1,4124):
Value_Matlab.append(Matlab_result[i+1][j])
X_Matlab.append(Matlab_mesh[j])
Phi_u0_Matlab = np.asarray(Value_Matlab)
plt.plot(X_Matlab,Phi_u0_Matlab, label = "Matlab")
### plot Python Results ###
Value_Python =[]
X_Python = []
tol = .001 # avoid hitting points outside the domain
x = np.linspace(0 + tol, 4 - tol, 20000)
points = [(x_, 0.01) for x_ in x]
Phi_u0_x1 = np.array([Phi_u0(point) for point in points])
plt.plot(x, Phi_u0_x1, label = "Python")
### Plot ###
plt.axis([0,4,-0.05,0.8])
plt.ylabel('Phi_u0')
plt.xlabel('x')
plt.legend()
plt.show()
cFile.close()
Solve()
```
%% Output
---------------------------------------------------------------------------
OSError Traceback (most recent call last)
<ipython-input-1-9bc1d91bd2c9> in <module>
228
229
--> 230 Solve()
<ipython-input-1-9bc1d91bd2c9> in Solve()
47 ### Import Matlab Results ###
48
---> 49 Matlab_result = np.loadtxt('resultsmatlab.txt')
50 Matlab_mesh = np.loadtxt('meshematlab.txt')
51
~/anaconda3/envs/fe38/lib/python3.8/site-packages/numpy/lib/npyio.py in loadtxt(fname, dtype, comments, delimiter, converters, skiprows, usecols, unpack, ndmin, encoding, max_rows)
979 fname = os_fspath(fname)
980 if _is_string_like(fname):
--> 981 fh = np.lib._datasource.open(fname, 'rt', encoding=encoding)
982 fencoding = getattr(fh, 'encoding', 'latin1')
983 fh = iter(fh)
~/anaconda3/envs/fe38/lib/python3.8/site-packages/numpy/lib/_datasource.py in open(path, mode, destpath, encoding, newline)
267
268 ds = DataSource(destpath)
--> 269 return ds.open(path, mode, encoding=encoding, newline=newline)
270
271
~/anaconda3/envs/fe38/lib/python3.8/site-packages/numpy/lib/_datasource.py in open(self, path, mode, encoding, newline)
621 encoding=encoding, newline=newline)
622 else:
--> 623 raise IOError("%s not found." % path)
624
625
OSError: resultsmatlab.txt not found.
%% Cell type:code id: tags:
``` python
```
......
Meshes/Like1DMesh.geo
+ 9
− 8
View file @ aef0d678
......@@ -3,19 +3,20 @@
lc = 10;
//+
Point(1) = {0, 0, 0, 0.1};
Point(1) = {0, 0, 0, lc/100};
//+
Point(2) = {0, 0.1, 0, 0.1};
Point(2) = {0, 0.01, 0, lc/100};
//+
Point(3) = {30, 0.1, 0, lc};
Point(3) = {30, 0.01, 0, lc};
//+
Point(4) = {30, 0, 0, lc};
//+
Point(5) = {1, 0.05, 0, 1};
//+
Point(6) = {1, 0, 0, 0.0001};
Point(6) = {1, 0, 0, 1};
......@@ -45,13 +46,13 @@ Physical Surface("Plane") = {1};
Field[1] = Distance;
//+
Field[1].NodesList = {6,5};
Field[1].NNodesByEdge = 1000;
Field[1].NNodesByEdge = 500;
Field[2] = Threshold;
Field[2].IField = 1;
Field[2].LcMin = lc /10000;
Field[2].LcMax = lc;
Field[2].LcMin = lc/1000;
Field[2].LcMax = lc/5;
Field[2].DistMin = 0.5;
Field[2].DistMax = 0.75;
Field[2].DistMax = 1;
Background Field = 2;
//+
......
comparison_FEM_FD.m
+ 11
− 4
View file @ aef0d678
%% Use below (this is now in radial cooordinates) for comparison with Python.
nu = 10^-15;
nu = 10^-6;
chi = 7/3;
b = nu^(1/3)*sqrt(chi/(chi-2));
e = sqrt(3/8*(chi-2));
T = Ternary_model(2, 'FRAP', [-1, b, 0.5, e, 1, 0.5], linspace(0.0, 400, 5000));
T.t = 0:0.1:20;
T.solve_tern_frap()
\ No newline at end of file
T = Ternary_model(0, 'FRAP', [-1, b, 0.5, e, 1, 0.5, 300], linspace(0.0, 400, 5000), 5);
T.t = 0:0.1:10;
T.solve_tern_frap()
%%
csvwrite('resultsmatlab.txt', T.sol);
csvwrite('meshematlab.txt', T.x);
csvwrite('phi_tot.txt', T.phi_t);
% writematrix(T.sol, 'resultsmatlab.txt');
% % csvwrite('meshematlab.txt', T.x);
\ No newline at end of file
like_1D.xmdf 0 → 100644
+ 1220
− 0
View file @ aef0d678
This diff is collapsed.
meshematlab.txt 0 → 100644
+ 1
− 0
View file @ aef0d678
This diff is collapsed.
phi_tot.txt 0 → 100644
+ 1
− 0
View file @ aef0d678
This diff is collapsed.
resultsmatlab.txt 0 → 100644
+ 101
− 0
View file @ aef0d678
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment