Commit 3796acbd authored by Jeremy BLEYER's avatar Jeremy BLEYER

Almost working nonlinear heat equation

parent afdaae79
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Feb 1 08:49:11 2019
@author: bleyerj
"""
from dolfin import *
import mfront_wrapper as mf
import numpy as np
length = 30e-3
width = 5.4e-3
mesh = RectangleMesh(Point(0., 0.), Point(length, width), 100, 10)
V = FunctionSpace(mesh, "CG", 1)
T = Function(V, name="Temperature")
def left(x, on_boundary):
return near(x[1], 0) and on_boundary
def right(x, on_boundary):
return near(x[0], length) and on_boundary
Tl = 300
Tr = 800
bc = [DirichletBC(V, Constant(Tl), left),
DirichletBC(V, Constant(Tr), right)]
facets = MeshFunction("size_t", mesh, 1)
ds = Measure("ds", subdomain_data=facets)
dt = Constant(1e-3)
theta = Constant(0.5)
material = mf.MFrontNonlinearMaterial("../materials/src/libBehaviour.so",
"HeatTransfer",
hypothesis="plane_strain")
problem = mf.MFrontNonlinearProblem(T, material, quadrature_degree=0)
problem.bc = bc
j = mf.ThermalFlux(T)
j.initialize_functions(mesh, 0)
problem.flux = j
H = mf.Flux(mf.Var(T), name="Enthalpy")
H.initialize_functions(mesh, 0)
problem.state_variables.append(H)
H0 = H.previous()
j0 = j.previous()
problem.L = (problem.u_*(H.function - H0)-dt*(theta*dot(grad(problem.u_), j.function)+
-(1-theta)*dot(grad(problem.u_), j0)))*problem.dx
problem.a = (problem.u_*H.tangents[0]*problem.du -
dt*theta*dot(grad(problem.u_), j.tangents[0]*grad(problem.du) +
j.tangents[1]*problem.du))*problem.dx
print(problem.material.data_manager.K.shape)
T.interpolate(Constant(Tl))
problem.solve(T.vector())
x = np.linspace(0, length)
import matplotlib.pyplot as plt
plt.plot(x, np.array([T(xi, width/2) for xi in x]))
#plt.figure()
#plot(project(flux.function[0], FunctionSpace(mesh, "DG",0)))
\ No newline at end of file
@DSL DefaultGenericBehaviour;
@Behaviour HeatTransfer;
@Author Thomas Helfer;
@Date 15 / 02 / 2019;
@Gradient TemperatureGradient gT;
gT.setGlossaryName("TemperatureGradient");
@Flux HeatFlux j;
j.setGlossaryName("HeatFlux");
@StateVariable real H;
H.setEntryName("Enthalpy");
@AdditionalTangentOperatorBlock dj_ddT;
@Parameter Tm = 933.15;
@Parameter ks = 210;
@Parameter Cs = 3.e6;
@Parameter kl = 95;
@Parameter Cl = 2.58e6;
@Parameter dHsl = 1.08048e9;
@LocalVariable thermalconductivity k;
// @LocalVariable real H0;
@LocalVariable real Ce;
@Integrator {
const auto T_ = T + dT;
// const auto Te = temperature(0.01);
// const auto Ts = Tm-Te;
// const auto Tl = Tm+Te;
// heat flux
k = (T_<Tm) ? ks : kl;
j = -k*(gT+dgT);
// enthalpy
if(T_<Tm){
Ce = Cs;
H = Cs*T_;
} else {
Ce = Cl;
H = Cl*(T_-Tm)+dHsl+Cs*Tm;
}
} // end of @Integrator
@TangentOperator {
dj_ddgT = -k * tmatrix<N, N, real>::Id();
dj_ddT = tvector<N, real>(real(0));
dH_ddT = Ce;
} // end of @TangentOperator
......@@ -6,9 +6,11 @@
@Gradient TemperatureGradient gT;
gT.setGlossaryName("TemperatureGradient");
@Flux HeatFlux j;
j.setGlossaryName("HeatFlux");
@AdditionalTangentOperatorBlock dj_ddT;
@Parameter A = 0.0375;
......
......@@ -173,12 +173,12 @@ private :
#line 17 "StationaryHeatTransfer.mfront"
#line 19 "StationaryHeatTransfer.mfront"
thermalconductivity k;
#line 14 "StationaryHeatTransfer.mfront"
#line 16 "StationaryHeatTransfer.mfront"
real A;
#line 15 "StationaryHeatTransfer.mfront"
#line 17 "StationaryHeatTransfer.mfront"
real B;
real minimal_time_step_scaling_factor;
real maximal_time_step_scaling_factor;
......@@ -315,11 +315,11 @@ using namespace tfel::math;
raise_if(smflag!=MechanicalBehaviour<MechanicalBehaviourBase::GENERALBEHAVIOUR,hypothesis,Type,false>::STANDARDTANGENTOPERATOR,
"invalid tangent operator flag");
bool computeTangentOperator_ = smt!=NOSTIFFNESSREQUESTED;
#line 20 "StationaryHeatTransfer.mfront"
#line 22 "StationaryHeatTransfer.mfront"
const auto T_ = this->T + this->dT;
#line 21 "StationaryHeatTransfer.mfront"
#line 23 "StationaryHeatTransfer.mfront"
this->k = 1 / (this->A + this->B * T_);
#line 22 "StationaryHeatTransfer.mfront"
#line 24 "StationaryHeatTransfer.mfront"
this->j = this->k * (this->gT + this->dgT);
this->updateIntegrationVariables();
this->updateStateVariables();
......@@ -364,9 +364,9 @@ bool computeConsistentTangentOperator(const SMType smt){
using namespace std;
using namespace tfel::math;
using std::vector;
#line 26 "StationaryHeatTransfer.mfront"
#line 28 "StationaryHeatTransfer.mfront"
this->dj_ddgT = this->k * tmatrix<N, N, real>::Id();
#line 27 "StationaryHeatTransfer.mfront"
#line 29 "StationaryHeatTransfer.mfront"
dj_ddT = -this->B * this->k * this->k * (this->gT + this->dgT);
return true;
}
......
......@@ -12,7 +12,7 @@ INCLUDES := -I../include \
CXXFLAGS := -Wall -Wfatal-errors -ansi $(shell tfel-config --oflags) -fPIC $(INCLUDES)
SRCCXX = IsotropicLinearHardeningPlasticity-generic.cxx IsotropicLinearHardeningPlasticity.cxx LogarithmicStrainPlasticity-generic.cxx LogarithmicStrainPlasticity.cxx StationaryHeatTransfer-generic.cxx StationaryHeatTransfer.cxx
SRCCXX = HeatTransfer-generic.cxx HeatTransfer.cxx IsotropicLinearHardeningPlasticity-generic.cxx IsotropicLinearHardeningPlasticity.cxx LogarithmicStrainPlasticity-generic.cxx LogarithmicStrainPlasticity.cxx MultiPhase-generic.cxx MultiPhase.cxx SHT3-generic.cxx SHT3.cxx StationaryHeatTransfer-generic.cxx StationaryHeatTransfer.cxx
makefiles1 = $(SRCCXX:.cxx=.d)
makefiles2 = $(makefiles1:.cpp=.d)
......@@ -22,7 +22,7 @@ makefiles = $(makefiles2)
all : libBehaviour.so
libBehaviour.so : IsotropicLinearHardeningPlasticity-generic.o IsotropicLinearHardeningPlasticity.o LogarithmicStrainPlasticity-generic.o LogarithmicStrainPlasticity.o StationaryHeatTransfer-generic.o StationaryHeatTransfer.o
libBehaviour.so : IsotropicLinearHardeningPlasticity-generic.o IsotropicLinearHardeningPlasticity.o LogarithmicStrainPlasticity-generic.o LogarithmicStrainPlasticity.o StationaryHeatTransfer-generic.o StationaryHeatTransfer.o HeatTransfer-generic.o HeatTransfer.o SHT3-generic.o SHT3.o MultiPhase-generic.o MultiPhase.o
@$(CXX) -shared $^ -o $@ -L"$(strip $(shell tfel-config --library-path))" $(patsubst %,-l%,$(shell tfel-config --library-dependency --material --mfront-profiling --physical-constants))
install :
......
......@@ -11,7 +11,13 @@ sources : {
"LogarithmicStrainPlasticity-generic.cxx",
"LogarithmicStrainPlasticity.cxx",
"StationaryHeatTransfer-generic.cxx",
"StationaryHeatTransfer.cxx"
"StationaryHeatTransfer.cxx",
"HeatTransfer-generic.cxx",
"HeatTransfer.cxx",
"SHT3-generic.cxx",
"SHT3.cxx",
"MultiPhase-generic.cxx",
"MultiPhase.cxx"
};
cppflags : {
"$(shell tfel-config --cppflags --compiler-flags)"
......@@ -40,7 +46,22 @@ epts : {
"StationaryHeatTransfer_Axisymmetrical",
"StationaryHeatTransfer_PlaneStrain",
"StationaryHeatTransfer_GeneralisedPlaneStrain",
"StationaryHeatTransfer_Tridimensional"
"StationaryHeatTransfer_Tridimensional",
"HeatTransfer_AxisymmetricalGeneralisedPlaneStrain",
"HeatTransfer_Axisymmetrical",
"HeatTransfer_PlaneStrain",
"HeatTransfer_GeneralisedPlaneStrain",
"HeatTransfer_Tridimensional",
"SHT3_AxisymmetricalGeneralisedPlaneStrain",
"SHT3_Axisymmetrical",
"SHT3_PlaneStrain",
"SHT3_GeneralisedPlaneStrain",
"SHT3_Tridimensional",
"MultiPhase_AxisymmetricalGeneralisedPlaneStrain",
"MultiPhase_Axisymmetrical",
"MultiPhase_PlaneStrain",
"MultiPhase_GeneralisedPlaneStrain",
"MultiPhase_Tridimensional"
};
};
headers : {
......@@ -55,6 +76,18 @@ headers : {
"MFront/GenericBehaviour/StationaryHeatTransfer-generic.hxx",
"TFEL/Material/StationaryHeatTransfer.hxx",
"TFEL/Material/StationaryHeatTransferBehaviourData.hxx",
"TFEL/Material/StationaryHeatTransferIntegrationData.hxx"
"TFEL/Material/StationaryHeatTransferIntegrationData.hxx",
"MFront/GenericBehaviour/HeatTransfer-generic.hxx",
"TFEL/Material/HeatTransfer.hxx",
"TFEL/Material/HeatTransferBehaviourData.hxx",
"TFEL/Material/HeatTransferIntegrationData.hxx",
"MFront/GenericBehaviour/SHT3-generic.hxx",
"TFEL/Material/SHT3.hxx",
"TFEL/Material/SHT3BehaviourData.hxx",
"TFEL/Material/SHT3IntegrationData.hxx",
"MFront/GenericBehaviour/MultiPhase-generic.hxx",
"TFEL/Material/MultiPhase.hxx",
"TFEL/Material/MultiPhaseBehaviourData.hxx",
"TFEL/Material/MultiPhaseIntegrationData.hxx"
};
};
......@@ -105,7 +105,13 @@ class Flux:
def update(self, x):
self.function.vector().set_local(x)
def previous(self):
try:
self.previous = Function(self.function.function_space(), name=self.name+"_previous")
self.previous.assign(self.function)
return self.previous
except:
raise ValueError("Function must be initialilzed first.")
class Stress(Flux):
def __init__(self, variable, name="sigma"):
......
......@@ -149,15 +149,28 @@ class MFrontNonlinearProblem(NonlinearProblem):
flattened_shapes = [s[0]*s[1] if len(s)==2 else s[0] for s in shapes]
for (i,t) in enumerate(self.flux.tangents):
s = flattened_shapes[i]
print(s, buff, K.shape, K[137,:])
# print(s, buff, K.shape, K[137,:])
t.vector().set_local(K[:,buff:buff+s].flatten())
#t.vector().set_local(K.flatten())
buff += s
print(self.flux.tangents[0].vector().get_local())
# sizes = self.material.get_state_variable_sizes()
# for (s, i) in self.state_variables:
# size = sizes[i]
# s.vector().set_local(self.material.data_manager.s1.internal_state_variables[:, i:(i+size)].flatten())
# print(self.flux.tangents[0].vector().get_local())
sizes = self.material.get_state_variable_sizes()
for i in range(len(self.state_variables)):
size = sizes[i]
s_values = self.material.data_manager.s1.internal_state_variables[:, i:(i+size)].flatten()
state = self.state_variables[i]
if isinstance(state, Flux):
state.function.vector().set_local(s_values)
shapes = [ufl.shape(t) for t in state.tangents]
print(shapes)
flattened_shapes = [s[0]*s[1] if len(s)==2 else s[0] if len(s)==1 else 1 for s in shapes]
# FIXME: Assumes gradients of state variables are defined after fluxes and in right order
for (j,t) in enumerate(state.tangents):
s = flattened_shapes[j]
t.vector().set_local(K[:,buff:buff+s].flatten())
buff += s
else:
state.vector().set_local(s_values)
def get_state_variable(self, name=None, position=None):
if name is not None:
......
......@@ -86,7 +86,7 @@ def gradient(g):
return nonsymmetric_tensor_to_vector(g)
def get_quadrature_element(cell, degree, dim=0):
if dim in [0, 1, (), (0,), (1,)]:
if dim in [0, 1, (), (0,), (1,), (0, 0), (0, 1), (1, 0)]:
return FiniteElement("Quadrature", cell, degree=degree,
quad_scheme='default')
elif type(dim) == int:
......
Markdown is supported
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