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

Readding Gamma_0 in Python and Matlab, still works.

parent 85254894
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
.gitignore
+ 1
− 0
View file @ fda5e299
......@@ -8,3 +8,4 @@
*.xdmf
*.xml
*.msh
*.txt
@Ternary_model/solve_tern_frap.m
+ 2
− 2
View file @ fda5e299
......@@ -17,8 +17,8 @@ function [c, f ,s] = flory_hugg_pde(x, t, u, dudx, a, b, e, u0,...
pt = Ternary_model.phi_tot(x, a, b, e, u0);
gra_a = Ternary_model.gradient_analytical(x, a, b, e);
% g0 = Ternary_model.gamma0(x, a, b, e_g0, u_g0);
g0 = 1;
g0 = Ternary_model.gamma0(x, a, b, e_g0, u_g0);
% g0 = 1;
c = 1;
f = g0*(1-pt)/pt*(pt*dudx-u*gra_a);
s = 0;
......
Like1D_Python_Matlab.ipynb
+ 11
− 5
View file @ fda5e299
%% 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 Solve():
set_log_level(40)
dt = 0.1
##PARAMETERS##
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 = 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('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)
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 = (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)
Form = (inner((Phi_u-Phi_u0)/dt, v) - inner(Gamma_0*((1-Phi_tot)/Phi_tot)*Phi_u*grad(Phi_tot) -
Gamma_0*(1-Phi_tot)*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 = XDMFFile('like_1D.xdmf')
cFile.write(Phi_u0, t)
ti = time.time()
### Computing ###
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)
cFile.close()
Solve()
```
%% Cell type:code id: tags:
``` python
```
......
comparison_FEM_FD.m
+ 1
− 1
View file @ fda5e299
......@@ -3,7 +3,7 @@ nu = 10^-6;
chi = 7/3;
b = nu^(1/3)*sqrt(chi/(chi-2));
e = sqrt(3/8*(chi-2));
T = Ternary_model(0, 'FRAP', [-1, b, 0.5, e, 1, 0.5, 300], linspace(0.0, 400, 5000), 5);
T = Ternary_model(0, 'FRAP', [-1, b, 0.5, e, 1, 0.5, 300], linspace(0.0, 400, 5000), 2);
T.t = 0:0.1:10;
T.solve_tern_frap()
%%
......
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