python-optipuls/run.py at develop · optipulsproject/python-optipuls · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
import argparse

import dolfin
from dolfin import Constant, as_matrix
import numpy as np

import optipuls.visualization as vis
from optipuls.simulation import Simulation
from optipuls.problem import Problem
import optipuls.material as material
import optipuls.optimization as optimization
from optipuls.time import TimeDomain
from optipuls.space import SpaceDomain
from optipuls.material import Material


# parse command line arguments
parser = argparse.ArgumentParser()
parser.add_argument('-o', '--output', default='../output')
parser.add_argument('-s', '--scratch', default='/scratch/OptiPuls/current')
args = parser.parse_args()


# set dolfin parameters
dolfin.set_log_level(40)
dolfin.parameters["form_compiler"]["quadrature_degree"] = 1


# set up the problem
problem = Problem()

time_domain = TimeDomain(0.020, 200)
problem.time_domain = time_domain


# experimental

mesh = dolfin.Mesh()
with dolfin.XDMFFile("examples/mesh/singlespot-XYZ/singlespot_XYZ.xdmf") as infile:
    infile.read(mesh)

mesh_boundaries = dolfin.MeshValueCollection("size_t", mesh, 2)
with dolfin.XDMFFile("examples/mesh/singlespot-XYZ/singlespot_XYZ_boundaries.xdmf") as infile:
    infile.read(mesh_boundaries, "ids")

space_domain = SpaceDomain(0.0050, 0.0002, 0.0005)
space_domain.dim = 3
space_domain._x = dolfin.SpatialCoordinate(mesh)
subdomain_data = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mesh_boundaries)
space_domain._ds = dolfin.Measure('ds', domain=mesh, subdomain_data=subdomain_data)

problem.space_domain = space_domain

#####

P_YAG = 3000.
absorb = 0.135
laser_pd = (absorb * P_YAG) / (np.pi * space_domain.R_laser**2)

problem.P_YAG = P_YAG
problem.laser_pd = laser_pd

problem.temp_amb = 295.
problem.implicitness = 1.
problem.convection_coeff = 20.
problem.radiation_coeff = 2.26 * 10**-9
problem.liquidus = 923.0
problem.solidus = 858.0

# optimization parameters
problem.beta_control = 10**2
problem.beta_velocity = 10**18
problem.velocity_max = 0.15
problem.beta_liquidity = 10**12
problem.beta_welding = 10**-2
problem.threshold_temp = 1000.
problem.target_point = dolfin.Point(0, 0, -.7 * space_domain.Z)
problem.pow_ = 20


# initialize FEM spaces
problem.V = dolfin.FunctionSpace(mesh, "CG", 1)
problem.V1 = dolfin.FunctionSpace(mesh, "DG", 0)

problem.theta_init = dolfin.project(problem.temp_amb, problem.V)

########################################


problem.material = Material.load('EN_AW-6082_T6.json')


# R = 0.0025;
# r = 0.0002;
# overlap = 0.6;
# d = 2 * r * (1 - overlap);


from optipuls.utils.laser import linear_rampdown
from optipuls.utils.io import save_as_pvd

print('Creating an initial simulation.')
control = np.zeros(time_domain.Nt)
simulation = Simulation(problem, control)

descent = optimization.gradient_descent(
        simulation, iter_max=20, step_init=2**-25)

vis.control_plot(
        descent[-1].control,
        labels=['Optimal Control'],
        outfile=args.scratch+'/optimal_control.png')
print(f'Gradient descent complete. See {args.scratch}/optimal_control.png')

# optimized = descent[-1]
# save_as_pvd(optimized.evo, problem.V, 'output/optimized.pvd')


# control = linear_rampdown(time_domain.timeline, t1=0.005, t2=0.010)
# simulation = Simulation(problem, control)
# print(simulation.evo.max())
# save_as_pvd(simulation.evo, problem.V, 'output/evo.pvd')
# theta_final = dolfin.Function(problem.V)
# theta_final.set_allow_extrapolation(True)
# theta_final.vector().set_local(simulation.evo[-1])

# from optipuls.utils.shiftedexpression import ShiftedExpression

# theta_shifted = ShiftedExpression(
#     expr_interpolate=(lambda p: theta_final(p)),
#     expr_extrapolate=(lambda _: problem.temp_amb),
#     shift=(d, 0, 0),
#     bounding_box_tree=mesh.bounding_box_tree(),
#     )

# theta_init_new = dolfin.Function(problem.V)
# theta_init_new = dolfin.interpolate(theta_shifted, problem.V)


# print('Creating a test simulation.')
# test_control = 0.5 + 0.1 * np.sin(0.5 * time_domain.timeline / np.pi)
# test_simulation = Simulation(problem, test_control)

# epsilons, deltas_fwd = optimization.gradient_test(
#         test_simulation, eps_init=10**-5, iter_max=15)
# vis.gradient_test_plot(
#         epsilons, deltas_fwd, outfile=args.scratch+'/gradient_test.png')
# print(f'Gradient test complete. See {args.scratch}/gradient_test.png')


# print('A linear rampdown simulation.')
# control = linear_rampdown(time_domain.timeline, t1=0.005, t2=0.010)
# simulation = Simulation(problem, control)
# save_as_pvd(simulation.evo, problem.V, 'output/evo.pvd')

# print('Creating an initial simulation.')
# control = np.zeros(time_domain.Nt)
# simulation = Simulation(problem, control)

# descent = optimization.gradient_descent(
#         simulation, iter_max=20, step_init=2**-25)

# vis.control_plot(
#         descent[-1].control,
#         labels=['Optimal Control'],
#         outfile=args.scratch+'/optimal_control.png')
# print(f'Gradient descent complete. See {args.scratch}/optimal_control.png')

# optimized = descent[-1]
# save_as_pvd(optimized.evo, problem.V, 'output/optimized.pvd')