Skip to content

Commit

Permalink
In work
Browse files Browse the repository at this point in the history
  • Loading branch information
@niravshah241
niravshah241 committed Dec 20, 2023
1 parent 9cb287d commit 2ba6236
Showing 1 changed file with 10 additions and 223 deletions.
233 changes: 10 additions & 223 deletions demo/thermomechanical_dlrbnicsx/dlrbnicsx_thermomechanical.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
import rbnicsx.online

from dlrbnicsx_thermomechanical_geometric_deformation import MeshDeformationWrapperClass
from dlrbnicsx_thermal import ThermalProblemOnDeformedDomain

from mpi4py import MPI
from petsc4py import PETSc
Expand All @@ -27,171 +28,6 @@
from dlrbnicsx.train_validate_test.train_validate_test import \
train_nn, validate_nn, online_nn, error_analysis

class ThermalProblemOnDeformedDomain(abc.ABC):
def __init__(self, mesh, subdomains, boundaries):
self._mesh = mesh
self._subdomains = subdomains
self._boundaries = boundaries
dx = ufl.Measure("dx", domain=mesh, subdomain_data=subdomains)
ds = ufl.Measure("ds", domain=self._mesh, subdomain_data=self._boundaries)
self.dx = dx
self._ds_sf = ds(11) + ds(20) + ds(21) + ds(22) + ds(23)
self._ds_bottom = ds(1) + ds(31)
self._ds_out = ds(30)
self._ds_sym = ds(5) + ds(9) + ds(12)
self._ds_top = ds(18) + ds(19) + ds(27) + ds(28) + ds(29)
x = ufl.SpatialCoordinate(self._mesh)
self._VT = dolfinx.fem.FunctionSpace(self._mesh, ("CG", 1))
uT, vT = ufl.TrialFunction(self._VT), ufl.TestFunction(self._VT)
self._trial, self._test = uT, vT
self._inner_product = ufl.inner(uT, vT) * x[0] * ufl.dx + \
ufl.inner(ufl.grad(uT), ufl.grad(vT)) * x[0] * ufl.dx
self.inner_product_action = \
rbnicsx.backends.bilinear_form_action(self._inner_product,
part="real")
self._T_f, self._T_out, self._T_bottom = 1773., 300., 300.
self._h_cf, self._h_cout, self._h_cbottom = 2000., 200., 200.
self._q_source = dolfinx.fem.Function(self._VT)
self._q_source.x.array[:] = 0.
self._q_top = dolfinx.fem.Function(self._VT)
self._q_top.x.array[:] = 0.
self.uT_func = dolfinx.fem.Function(self._VT)
self.mu_ref = [0.6438, 0.4313, 1., 0.5]

self._max_iterations = 20
self.rtol = 1.e-4
self.atol = 1.e-12

def thermal_diffusivity_1(self, sym_T):
conditions = [ufl.le(sym_T, 293.), ufl.And(ufl.ge(sym_T, 293.), ufl.le(sym_T, 673.)), ufl.And(ufl.ge(sym_T, 673.), ufl.le(sym_T, 1273.)), ufl.ge(sym_T, 1273.)]
interps = [8.17484662576687e-6*sym_T**2 - 0.00926193251533741*sym_T + 18.0819438190184, 8.17484662576687e-6*sym_T**2 - 0.00926193251533741*sym_T + 18.0819438190184, 1.76073619631904e-6*sym_T**2 - 0.000628539877300632*sym_T + 15.176807196319, 1.76073619631904e-6*sym_T**2 - 0.000628539877300632*sym_T + 15.176807196319]
assert len(conditions) == len(interps)
d_func = ufl.conditional(conditions[0], interps[0], interps[0])
for i in range(1, len(conditions)):
d_func = ufl.conditional(conditions[i], interps[i], d_func)
return d_func

def thermal_diffusivity_2(self, sym_T):
conditions = [ufl.le(sym_T, 293.), ufl.And(ufl.ge(sym_T, 293.), ufl.le(sym_T, 693.)), ufl.And(ufl.ge(sym_T, 673.), ufl.le(sym_T, 1273.)), ufl.ge(sym_T, 1273.)]
interps = [0.000299054192229039*sym_T**2 - 0.36574217791411*sym_T + 130.838954780164, 0.000299054192229039*sym_T**2 - 0.36574217791411*sym_T + 130.838954780164, -1.10434560327197e-5*sym_T**2 + 0.0516492566462166*sym_T - 9.61326294938646, -1.10434560327197e-5*sym_T**2 + 0.0516492566462166*sym_T - 9.61326294938646]
assert len(conditions) == len(interps)
d_func = ufl.conditional(conditions[0], interps[0], interps[0])
for i in range(1, len(conditions)):
d_func = ufl.conditional(conditions[i], interps[i], d_func)
return d_func

def thermal_diffusivity_3(self, sym_T):
conditions = [ufl.le(sym_T, 293.), ufl.And(ufl.ge(sym_T, 293.), ufl.le(sym_T, 693.)), ufl.And(ufl.ge(sym_T, 673.), ufl.le(sym_T, 1273.)), ufl.ge(sym_T, 1273.)]
interps = [8.17484662576687e-6*sym_T**2 - 0.00926193251533741*sym_T + 18.0819438190184, 8.17484662576687e-6*sym_T**2 - 0.00926193251533741*sym_T + 18.0819438190184, 1.76073619631904e-6*sym_T**2 - 0.000628539877300632*sym_T + 15.176807196319, 1.76073619631904e-6*sym_T**2 - 0.000628539877300632*sym_T + 15.176807196319]
assert len(conditions) == len(interps)
d_func = ufl.conditional(conditions[0], interps[0], interps[0])
for i in range(1, len(conditions)):
d_func = ufl.conditional(conditions[i], interps[i], d_func)
return d_func

def thermal_diffusivity_4(self, sym_T):
conditions = [ufl.le(sym_T, 293.), ufl.And(ufl.ge(sym_T, 293.), ufl.le(sym_T, 693.)), ufl.And(ufl.ge(sym_T, 673.), ufl.le(sym_T, 1273.)), ufl.ge(sym_T, 1273.)]
interps = [8.17484662576687e-6*sym_T**2 - 0.00926193251533741*sym_T + 18.0819438190184, 8.17484662576687e-6*sym_T**2 - 0.00926193251533741*sym_T + 18.0819438190184, 1.76073619631904e-6*sym_T**2 - 0.000628539877300632*sym_T + 15.176807196319, 1.76073619631904e-6*sym_T**2 - 0.000628539877300632*sym_T + 15.176807196319]
assert len(conditions) == len(interps)
d_func = ufl.conditional(conditions[0], interps[0], interps[0])
for i in range(1, len(conditions)):
d_func = ufl.conditional(conditions[i], interps[i], d_func)
return d_func

def thermal_diffusivity_5(self, sym_T):
conditions = [ufl.le(sym_T, 293.), ufl.And(ufl.ge(sym_T, 293.), ufl.le(sym_T, 693.)), ufl.And(ufl.ge(sym_T, 673.), ufl.le(sym_T, 1273.)), ufl.ge(sym_T, 1273.)]
interps = [3.08018064076346e-5*sym_T**2 - 0.0376497392638036*sym_T + 31.7270693260054, 3.08018064076346e-5*sym_T**2 - 0.0376497392638036*sym_T + 31.7270693260054, -2.79311520109062e-6*sym_T**2 + 0.00756902522154049*sym_T + 16.5109550766871, -2.79311520109062e-6*sym_T**2 + 0.00756902522154049*sym_T + 16.5109550766871]
assert len(conditions) == len(interps)
d_func = ufl.conditional(conditions[0], interps[0], interps[0])
for i in range(1, len(conditions)):
d_func = ufl.conditional(conditions[i], interps[i], d_func)
return d_func

def thermal_diffusivity_6(self, sym_T):
conditions = [ufl.le(sym_T, 293.), ufl.And(ufl.ge(sym_T, 293.), ufl.le(sym_T, 693.)), ufl.And(ufl.ge(sym_T, 673.), ufl.le(sym_T, 1273.)), ufl.ge(sym_T, 1273.)]
interps = [8.17484662576687e-6*sym_T**2 - 0.00926193251533741*sym_T + 18.0819438190184, 8.17484662576687e-6*sym_T**2 - 0.00926193251533741*sym_T + 18.0819438190184, 1.76073619631904e-6*sym_T**2 - 0.000628539877300632*sym_T + 15.176807196319, 1.76073619631904e-6*sym_T**2 - 0.000628539877300632*sym_T + 15.176807196319]
assert len(conditions) == len(interps)
d_func = ufl.conditional(conditions[0], interps[0], interps[0])
for i in range(1, len(conditions)):
d_func = ufl.conditional(conditions[i], interps[i], d_func)
return d_func

def thermal_diffusivity_7(self, sym_T):
conditions = [ufl.le(sym_T, 293.), ufl.And(ufl.ge(sym_T, 293.), ufl.le(sym_T, 693.)), ufl.And(ufl.ge(sym_T, 673.), ufl.le(sym_T, 1273.)), ufl.ge(sym_T, 1273.)]
interps = [0.000299054192229039*sym_T**2 - 0.36574217791411*sym_T + 130.838954780164, 0.000299054192229039*sym_T**2 - 0.36574217791411*sym_T + 130.838954780164, -1.10434560327197e-5*sym_T**2 + 0.0516492566462166*sym_T - 9.61326294938646, -1.10434560327197e-5*sym_T**2 + 0.0516492566462166*sym_T - 9.61326294938646]
assert len(conditions) == len(interps)
d_func = ufl.conditional(conditions[0], interps[0], interps[0])
for i in range(1, len(conditions)):
d_func = ufl.conditional(conditions[i], interps[i], d_func)
return d_func

@property
def lhs_form(self):
uT_func, vT = self.uT_func, self._test
x = ufl.SpatialCoordinate(self._mesh)
dx = self.dx
a_T = \
ufl.inner(self.thermal_diffusivity_1(uT_func) * ufl.grad(uT_func), ufl.grad(vT)) * x[0] * dx(1) + \
ufl.inner(self.thermal_diffusivity_2(uT_func) * ufl.grad(uT_func), ufl.grad(vT)) * x[0] * dx(2) + \
ufl.inner(5.3 * ufl.grad(uT_func), ufl.grad(vT)) * x[0] * dx(3) + \
ufl.inner(4.75 * ufl.grad(uT_func), ufl.grad(vT)) * x[0] * dx(4) + \
ufl.inner(self.thermal_diffusivity_5(uT_func) * ufl.grad(uT_func), ufl.grad(vT)) * x[0] * dx(5) + \
ufl.inner(45.6 * ufl.grad(uT_func), ufl.grad(vT)) * x[0] * dx(6) + \
ufl.inner(self.thermal_diffusivity_7(uT_func) * ufl.grad(uT_func), ufl.grad(vT)) * x[0] * dx(7) + \
ufl.inner(self._h_cf * uT_func, vT) * x[0] * self._ds_sf + \
ufl.inner(self._h_cout * uT_func, vT) * x[0] * self._ds_out + \
ufl.inner(self._h_cbottom * uT_func, vT) * x[0] * self._ds_bottom
return a_T

@property
def rhs_form(self):
vT = self._test
x = ufl.SpatialCoordinate(self._mesh)
dx = self.dx
l_T = \
ufl.inner(self._q_source, vT) * x[0] * dx + \
self._h_cf * vT * self._T_f * x[0] * self._ds_sf + \
self._h_cout * vT * self._T_out * x[0] * self._ds_out + \
self._h_cbottom * vT * self._T_bottom * x[0] * self._ds_bottom - \
ufl.inner(self._q_top, vT) * x[0] * self._ds_top
return l_T

@property
def set_problem(self):
problemNonlinear = \
NonlinearProblem(self.lhs_form - self.rhs_form,
self.uT_func, bcs=[])
return problemNonlinear

def solve(self, mu):
vT = self._test
self.mu = mu
self.uT_func.x.array[:] = 350.
self.uT_func.x.scatter_forward()
problemNonlinear = self.set_problem
solution = dolfinx.fem.Function(self._VT)
with MeshDeformationWrapperClass(self._mesh, self._boundaries,
self.mu_ref, self.mu):
solver = NewtonSolver(self._mesh.comm, problemNonlinear)
solver.convergence_criterion = "incremental"

solver.rtol = 1e-10
solver.report = True
ksp = solver.krylov_solver
ksp.setFromOptions()
# dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)

n, converged = solver.solve(self.uT_func)
# print(f"Computed solution array: {uT_func.x.array}")
assert (converged)
print(f"Number of interations: {n:d}")

solution.x.array[:] = self.uT_func.x.array.copy()
solution.x.scatter_forward()
x = ufl.SpatialCoordinate(self._mesh)
print(f"Temperature field norm: {self._mesh.comm.allreduce(dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(solution, solution) * x[0] * ufl.dx + ufl.inner(ufl.grad(solution), ufl.grad(solution)) * x[0] * ufl.dx)))}")
return solution

class MechanicalProblemOnDeformedDomain(abc.ABC):
def __init__(self, mesh, subdomains, boundaries, thermalproblem):
self._mesh = mesh
Expand Down Expand Up @@ -404,57 +240,6 @@ def solve(self, mu):
print(f"Displacement field norm: {self._mesh.comm.allreduce(dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(uM_func, uM_func) * x[0] * ufl.dx + ufl.inner(self.epsilon(uM_func), self.epsilon(uM_func)) * x[0] * ufl.dx)))}")
return uM_func

class ThermalPODANNReducedProblem(abc.ABC):
def __init__(self, thermal_problem) -> None:
self._basis_functions = rbnicsx.backends.FunctionsList(thermal_problem._VT)
uT, vT = ufl.TrialFunction(thermal_problem._VT), ufl.TestFunction(thermal_problem._VT)
x = ufl.SpatialCoordinate(thermal_problem._mesh)
self._inner_product = ufl.inner(uT, vT) * x[0] * ufl.dx + \
ufl.inner(ufl.grad(uT), ufl.grad(vT)) * x[0] * ufl.dx
self._inner_product_action = \
rbnicsx.backends.bilinear_form_action(self._inner_product,
part="real")
self.input_scaling_range = [-1., 1.]
self.output_scaling_range = [-1., 1.]
self.input_range = \
np.array([[0.55, 0.35, 0.8, 0.4], [0.75, 0.55, 1.2, 0.6]])
self.output_range = [-6., 3.]
self.regularisation = "EarlyStopping"

def reconstruct_solution(self, reduced_solution):
"""Reconstructed reduced solution on the high fidelity space."""
return self._basis_functions[:reduced_solution.size] * \
reduced_solution

def compute_norm(self, function):
"""Compute the norm of a function inner product
on the reference domain."""
return np.sqrt(self._inner_product_action(function)(function))

def project_snapshot(self, solution, N):
return self._project_snapshot(solution, N)

def _project_snapshot(self, solution, N):
projected_snapshot = rbnicsx.online.create_vector(N)
A = rbnicsx.backends.\
project_matrix(self._inner_product_action,
self._basis_functions[:N])
F = rbnicsx.backends.\
project_vector(self._inner_product_action(solution),
self._basis_functions[:N])
ksp = PETSc.KSP()
ksp.create(projected_snapshot.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.setFromOptions()
ksp.solve(F, projected_snapshot)
return projected_snapshot

def norm_error(self, u, v):
return self.compute_norm(u-v)/self.compute_norm(u)


class MechanicalPODANNReducedProblem(abc.ABC):
def __init__(self, mechanical_problem) -> None:
self._basis_functions = rbnicsx.backends.FunctionsList(mechanical_problem._VM)
Expand Down Expand Up @@ -525,6 +310,7 @@ def norm_error(self, u, v):
thermal_problem_parametric = \
ThermalProblemOnDeformedDomain(mesh, cell_tags, facet_tags)

'''
# solution_mu = thermal_problem_parametric.solve(mu_ref)
# print(f"Solution norm at mu:{mu_ref}: {thermal_problem_parametric.inner_product_action(solution_mu)(solution_mu)}")
Expand All @@ -540,7 +326,7 @@ def norm_error(self, u, v):
with dolfinx.io.XDMFFile(mesh.comm, computed_file, "w") as solution_file:
solution_file.write_mesh(mesh)
solution_file.write_function(uT_func_plot)

'''
# Thermal POD Starts ###

def generate_training_set(samples=pod_samples):
Expand All @@ -554,7 +340,7 @@ def generate_training_set(samples=pod_samples):
training_set_3)))
return training_set


'''
thermal_training_set = rbnicsx.io.on_rank_zero(mesh.comm, generate_training_set)
Nmax = 30
Expand Down Expand Up @@ -614,7 +400,7 @@ def generate_training_set(samples=pod_samples):
plt.savefig("thermal_eigenvalues.png")
print(f"Eigenvalues (Thermal): {thermal_positive_eigenvalues}")

'''
# Thermal POD Ends ###

# 5. ANN implementation
Expand Down Expand Up @@ -648,7 +434,7 @@ def generate_ann_output_set(problem, reduced_problem,
len(reduced_problem._basis_functions)).array.astype("f")
return output_set


'''
# Training dataset
thermal_ann_input_set = generate_ann_input_set(samples=ann_samples)
# np.random.shuffle(thermal_ann_input_set)
Expand Down Expand Up @@ -753,7 +539,7 @@ def generate_ann_output_set(problem, reduced_problem,
print(f"Error: {thermal_error_numpy[i]}")
# ### Thermal ANN ends

'''
mechanical_problem_parametric = \
MechanicalProblemOnDeformedDomain(mesh, cell_tags, facet_tags,
thermal_problem_parametric)
Expand Down Expand Up @@ -943,6 +729,7 @@ def generate_ann_output_set(problem, reduced_problem,

# ### Online phase ###
online_mu = np.array([0.45, 0.56, 0.9, 0.7])
'''
thermal_fem_solution = thermal_problem_parametric.solve(online_mu)
thermal_rb_solution = \
thermal_reduced_problem.reconstruct_solution(
Expand Down Expand Up @@ -998,7 +785,7 @@ def generate_ann_output_set(problem, reduced_problem,
print(thermal_reduced_problem.norm_error(thermal_fem_solution, thermal_rb_solution))
print(thermal_reduced_problem.compute_norm(thermal_error_function))

'''
mechanical_fem_solution = mechanical_problem_parametric.solve(online_mu)
mechanical_rb_solution = \
mechanical_reduced_problem.reconstruct_solution(
Expand Down Expand Up @@ -1051,5 +838,5 @@ def generate_ann_output_set(problem, reduced_problem,
print(mechanical_reduced_problem.norm_error(mechanical_fem_solution, mechanical_rb_solution))
print(mechanical_reduced_problem.compute_norm(mechanical_error_function))

print(f"Training time (Thermal): {thermal_elapsed_time}")
# print(f"Training time (Thermal): {thermal_elapsed_time}")
print(f"Training time (Mechanical): {mechanical_elapsed_time}")

0 comments on commit 2ba6236

Please sign in to comment.