Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

update #1357

Merged
merged 1 commit into from
Nov 21, 2024
Merged

update #1357

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
112 changes: 52 additions & 60 deletions app/tssim/ocp_opt/jovan/ocp_opt.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,12 +15,12 @@
logger.setLevel('ERROR') #积分子问题

bm.set_backend("numpy")
pde = example_1()
n = 20
q = 4
pde = example_1(c=0)
n = 10
q = 3
T = 1
nt = 30
maxit = 3
nt = 200
maxit = 15

mesh = TriangleMesh.from_box(pde.domain(), nx=n, ny=n)
timeline = UniformTimeLine(0, T, nt)
Expand All @@ -29,63 +29,53 @@
yspace= LagrangeFESpace(mesh, p=1)
space = LagrangeFESpace(mesh, p=1)
pspace = TensorFunctionSpace(space, (2,-1))
solver = ocp_opt_solver(mesh, yspace, pspace, pde, timeline)
solver = ocp_opt_solver(mesh, yspace, pspace, pde, timeline, q=q)

ygodf = yspace.number_of_global_dofs()
ygdof = yspace.number_of_global_dofs()
pgdof = pspace.number_of_global_dofs()
yisbdof = yspace.is_boundary_dof()
pisbdof = pspace.is_boundary_dof()
isbdof = bm.concatenate([yisbdof, pisbdof], axis=0)#(362,)


ally = [None]*(nt+1)
allp = [None]*(nt+1)
allu = [None]*(nt+1)
allz = [None]*(nt+1)

y0 = yspace.function(yspace.interpolate(partial(pde.y_solution, time=0)))
y0t = yspace.interpolate(partial(pde.y_t_solution, time=0))
p0 = pspace.function(pspace.interpolate(partial(pde.p_solution, time=0))) #(242,) p0x1,p0x2
p0 = pspace.function(pspace.interpolate(partial(pde.p_solution, time=0)))
zn = yspace.function(yspace.interpolate(partial(pde.z_solution, time=T)))
znt = yspace.interpolate(partial(pde.z_t_solution, time=0))

zn0 = yspace.function()
zn0[:] = zn[:]
zn1 = yspace.function()
zn2 = yspace.function()
q2 = pspace.function()

un = yspace.function()
un[:] = solver.z_to_u(zn) #积分子
allu[-1] = un
ally[0] = y0
allp[0] = p0
allz[-1] = zn

A0 = solver.Forward_BForm_A0().assembly()
b0_LForm = solver.Forward_LForm_b0()

A = solver.Forward_BForm_A().assembly()
FA = solver.Forward_BForm_A().assembly()
Forward_b_LForm = solver.Forward_LForm_b()

An = solver.Forward_BForm_A0().assembly()
bn_LForm = solver.Backward_LForm_bn()

Backward_A = solver.Forward_BForm_A().assembly()
BA = solver.Forward_BForm_A().assembly()
Backward_b_LForm = solver.Backward_LForm_b()

for k in range(maxit):
y1 = yspace.function()
p1 = pspace.function()


## 正向求解第0步
solver.Forward_0_update(allu[1])
b0 = b0_LForm.assembly()

Fb0 = b0_LForm.assembly()
BC = DirichletBC(space=(yspace,pspace), threshold=(None, None),
gd=(partial(pde.y_solution,time=dt),
partial(pde.p_solution,time=dt)),method='interp')
A0, b0 = BC.apply(A0, b0)
x0 = spsolve(A0, b0, solver='scipy')
A0, b0 = BC.apply(A0, Fb0)
x0 = spsolve(A0, b0, solver='mumps')

y1[:] = x0[:ygodf]
y1[:] = x0[:ygdof]
p1[:] = x0[-pgdof:]
ally[1] = y1
allp[1] = p1
Expand All @@ -101,15 +91,15 @@
p2 = pspace.function()

solver.Forward_update(ally[t1index-1], ally[t1index], allu[t1index+1], t1+dt)
Forward_b = Forward_b_LForm.assembly()
Fb = Forward_b_LForm.assembly()

BC = DirichletBC(space=(yspace,pspace), threshold=(None, None),
gd=(partial(pde.y_solution,time=t1+dt),
partial(pde.p_solution,time=t1+dt)),method='interp')
Forward_A, Forward_b = BC.apply(A, Forward_b)
Forward_x = spsolve(Forward_A, Forward_b, solver='scipy')
Forward_A, Forward_b = BC.apply(FA, Fb)
Forward_x = spsolve(Forward_A, Forward_b, solver='mumps')

y2[:] = Forward_x[:ygodf]
y2[:] = Forward_x[:ygdof]
p2[:] = Forward_x[-pgdof:]
ally[t1index+1] = y2
allp[t1index+1] = p2
Expand All @@ -119,52 +109,54 @@
t1 = timeline.current_time_level()
t1index = timeline.current_time_level_index()

zn1 = yspace.function()
solver.Backward_n_update(ally[t1index-1], allp[t1index-1])
bn = bn_LForm.assembly()

BC = DirichletBC(space=(yspace,pspace), threshold=(None, None),
gd=(partial(pde.y_solution,time = T-dt),
partial(pde.p_solution,time = T-dt)),method='interp')
gd=(partial(pde.z_solution,time = T-dt),
partial(pde.q_solution,time = T-dt)),method='interp')
An, bn = BC.apply(An, bn)
xn = spsolve(An, bn, solver='scipy')
zn1[:] = xn[:ygodf]
xn = spsolve(An, bn, solver='mumps')
zn1[:] = xn[:ygdof]
q2[:] = xn[-pgdof:]

un1 = yspace.function()
un1[:] = solver.z_to_u(zn1)
allu[t1index-1] = un1
allz[t1index-1] = zn1

timeline.backward()

## 反向求解
for i in bm.arange(nt-1, 0, -1):
t1 = timeline.current_time_level()
t1index = timeline.current_time_level_index()

zn2 = yspace.function()

##求第i-1步的z,q
solver.Backward_update(zn0, zn1, ally[t1index-1], allp[t1index-1], t1-dt)
Backward_b = Backward_b_LForm.assembly()
solver.Backward_update(allz[t1index+1], allz[t1index], ally[t1index-1], allp[t1index-1], t1-dt)
Bb = Backward_b_LForm.assembly()
BC = DirichletBC(space=(yspace,pspace), threshold=(None, None),
gd=(partial(pde.y_solution,time = t1-dt),
partial(pde.p_solution,time = t1-dt)),method='interp')
Backward_A, Backward_b = BC.apply(Backward_A, Backward_b)
Backward_x = spsolve(Backward_A, Backward_b, solver='scipy')
gd=(partial(pde.z_solution,time = t1-dt),
partial(pde.q_solution,time = t1-dt)),method='interp')
Backward_A, Backward_b = BC.apply(BA, Bb)
Backward_x = spsolve(Backward_A, Backward_b, solver='mumps')

zn2[:] = Backward_x[:ygodf]
zn2[:] = Backward_x[:ygdof]
q2 = Backward_x[-pgdof:]

##求第i-1步的u
un2 = yspace.function()
un2[:] = solver.z_to_u(zn2)
allu[t1index-1] = un2

zn0[:] = zn1[:]
zn1[:] = zn2[:]
allz[t1index-1] = zn2
timeline.backward()

ysolution = yspace.function(yspace.interpolate(partial(pde.y_solution, time=T)))
z_bar = solver.solve_z_bar(allz)
un_pre = allu[0]
if un_pre == None:
un_pre = 0
for i in range(nt+1):
ufunction = yspace.function()
ufunction[:] = bm.max((0,z_bar)) - allz[i]
allu[i] = ufunction
print(f"第{k}次的前后u差别",bm.sum(allu[0] - un_pre))

#ysolution = yspace.function(yspace.interpolate(partial(pde.y_solution, time=T)))
usolution = yspace.function(yspace.interpolate(partial(pde.u_solution, time=0)))
errory = mesh.error(ally[-1], ysolution)
#errory = mesh.error(ally[-1], ysolution)
erroru = mesh.error(allu[0], usolution)
print(errory)
#print(errory)
print(erroru)
19 changes: 7 additions & 12 deletions app/tssim/ocp_opt/jovan/solver_update.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@
from fealpy.backend import backend_manager as bm
from fealpy.sparse import COOTensor
from functools import partial
from tensor_mass_integrator import TensorMassIntegrator
from ocp_opt_pde import example_1

class ocp_opt_solver():
Expand Down Expand Up @@ -89,7 +88,6 @@ def coef(p, index=None):
Ly.add_integrator(self.forward_0_b_coef)

Lp = LinearForm(pspace)
#TODO:检查是不是一直是0
L = LinearBlockForm([Ly, Lp])
return L

Expand All @@ -103,7 +101,7 @@ def coef_u1(bcs, index=None):
result = (dt**2)*u1(bcs)
return result
self.forward_0_b_coef.source = coef_u1

self.forward_0_b_coef.clear()

def Forward_LForm_b(self):
yspace = self.yspace
Expand Down Expand Up @@ -145,12 +143,6 @@ def coef_c(p, index=None):
self.c_coef.clear()


### 反向求解
def z_to_u(self, z1):
result = bm.max(self.mesh.integral(z1), 0) - z1 #积分子
return result

## TODO:Pd Yd 的求解
def Backward_LForm_bn(self):
yspace = self.yspace
pspace = self.pspace
Expand Down Expand Up @@ -253,6 +245,9 @@ def coef_b_q(bcs, index=None):
self.backward_q_b_coef.source = coef_b_q
self.backward_q_b_coef.clear()




def solve_z_bar(self, allz):
dt = self.dt
integral_z = bm.array([self.mesh.integral(i, q=self.q) for i in allz],dtype=bm.float64)
z_bar = (dt/2)*(integral_z[:-1] + integral_z[1:])
return bm.sum(z_bar)

Loading