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

Jump to bottom

modify the device option in cmsmooth fem #1236

Merged
merged 3 commits into from
Oct 18, 2024
Merged

modify the device option in cmsmooth fem #1236

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
4e65665
functioninit
Oct 18, 2024
08d12fe
Merge remote-tracking branch 'upstream/master'
Oct 18, 2024
800cc29
光滑元算例
Oct 18, 2024
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
20 changes: 12 additions & 8 deletions example/experiment/fem/doublelaplace_fem_dirichlet.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
from fealpy import logger
from fealpy.solver import spsolve
logger.setLevel('INFO')
## 参数解析
parser = argparse.ArgumentParser(description=
Expand All @@ -34,7 +35,7 @@
help='初始网格剖分段数.')

parser.add_argument('--maxit',
default=2, type=int,
default=4, type=int,
help='默认网格加密求解的次数, 默认加密求解 4 次')

parser.add_argument('--backend',
Expand All @@ -60,8 +61,8 @@
x = sp.symbols('x')
y = sp.symbols('y')
u = (sp.sin(2*sp.pi*y)*sp.sin(2*sp.pi*x))**2
pde = DoubleLaplacePDE(u)
ulist = get_flist(u)[:3]
pde = DoubleLaplacePDE(u, device=device)
ulist = get_flist(u, device=device)[:3]
mesh = TriangleMesh.from_box([0,1,0,1], n, n, device=device)

ikwargs = bm.context(mesh.cell)
Expand Down Expand Up @@ -105,16 +106,19 @@
gdof = space.number_of_global_dofs()
NDof[i] = 1/4/2**i
bc1 = DirichletBC(space, gd = ulist)
#import ipdb
#ipdb.set_trace()
A, F = bc1.apply(A, F)
tmr.send(f'第{i}次边界处理时间')
A = A.to_scipy()
#A = A.to_scipy()

from numpy.linalg import cond
print(gdof)
print(cond(A.toarray()))
#from numpy.linalg import cond
#print(gdof)
#print(cond(A.toarray()))
#A = coo_matrix(A)
#A = csr_matrix((A.values(), A.indices()),A.shape)
uh[:] = bm.tensor(spsolve(A, F))
#uh[:] = bm.tensor(spsolve(A, F))
uh[:] = spsolve(A, F, "cupy")

#uh[:] = cg(A, F, maxiter=400000, atol=1e-14, rtol=1e-14)
tmr.send(f'第{i}次求解器时间')
Expand Down
2 changes: 1 addition & 1 deletion fealpy/functionspace/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,4 +13,4 @@
from .first_nedelec_fe_space_3d import FirstNedelecFiniteElementSpace3d

from.second_nedelec_fe_space_2d import SecondNedelecFiniteElementSpace2d
from.second_nedelec_fe_space_3d import SecondNedelecFiniteElementSpace3d
from.second_nedelec_fe_space_3d import SecondNedelecFiniteElementSpace3d
10 changes: 5 additions & 5 deletions fealpy/functionspace/cm_conforming_fe_space.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -343,7 +343,7 @@ def boundary_interpolate(self, gD, uh, threshold=None):
for r in range(1, 2*m+1):
val = gD[r](node)
multiIndex = self.mesh.multi_index_matrix(r, 1)
symidx, num = symmetry_index(2, r)
symidx, num = symmetry_index(2, r, device=self.device)

#idx = self.mesh.multi_index_matrix(r, 1)
#num = factorial(r)/bm.prod(factorial(idx), axis=1)
Expand Down Expand Up @@ -371,7 +371,7 @@ def boundary_interpolate(self, gD, uh, threshold=None):
bcoeff = bm.einsum('el, il->ei', ffval, b2l)
uh[e2id[:, k:l]] = bcoeff[:, 2*m+1-r: -2*m-1+r]
else:
symidx, num = symmetry_index(2, r)
symidx, num = symmetry_index(2, r, device=self.device)
nnn = span_array(n[:, None, :], bm.array([r]))
nnn = nnn.reshape(NE, -1)[:, symidx]

Expand Down Expand Up @@ -427,9 +427,9 @@ def interpolation(self, flist):
bcoeff = bm.einsum('el, il->ei', ffval, b2l)
fI[e2id[:, k:l]] = bcoeff[:, 2*m+1-r: -2*m-1+r]
else:
symidx, num = symmetry_index(2, r)
symidx, num = symmetry_index(2, r, device=self.device)

nnn = span_array(n[:, None, :], bm.array([r]))
nnn = span_array(n[:, None, :], bm.array([r], device=self.device))
nnn = nnn.reshape(NE, -1)[:, symidx]

#GD = self.mesh.geo_dimension()
Expand Down Expand Up @@ -480,7 +480,7 @@ def interpolation(self, flist):
for r in range(1, 2*m+1):
val = flist[r](node)
multiIndex = self.mesh.multi_index_matrix(r, 1)
symidx, num = symmetry_index(2, r)
symidx, num = symmetry_index(2, r, device=self.device)

#idx = self.mesh.multi_index_matrix(r, 1)
#num = factorial(r)/bm.prod(factorial(idx), axis=1)
Expand Down
2 changes: 1 addition & 1 deletion fealpy/functionspace/functional.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -144,7 +144,7 @@ def symmetry_index(d, r, dtype=None, device=None):

P = bm.concatenate([bm.tensor([1],device=device), bm.cumprod(bm.arange(r+1, device=device)[1:], axis=0)],
axis=0, dtype=dtype)
num = P[r]/bm.prod(P[midx], axis=1)
num = P[r]/bm.prod(P[midx], axis=1, dtype=bm.float64)
return symidx, num

def multi_index2d_to_index(midx):
Expand Down
77 changes: 48 additions & 29 deletions fealpy/pde/biharmonic_triharmonic_2d.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,23 +97,26 @@ def source(self, p):
#sin = bm.sin
#pi = bm.pi
#cos = bm.cos
print(bm.array(self.L2u(x_cpu, y_cpu),device=self.device))
return bm.array(self.L2u(x_cpu, y_cpu),device=self.device)

def solution(self, p):
x = p[..., 0]
y = p[..., 1]
return self.u(x, y)
x = bm.to_numpy(x)
y = bm.to_numpy(y)
return bm.array(self.u(x, y), device=self.device)

def gradient(self, p):
x = p[..., 0]
y = p[..., 1]
#sin = bm.sin
#pi = bm.pi
#cos = bm.cos
val = bm.zeros(p.shape, dtype=bm.float64)
val[..., 0] = self.ux(x, y)
val[..., 1] = self.uy(x, y)
x = bm.to_numpy(x)
y = bm.to_numpy(y)
val = bm.zeros(p.shape, dtype=bm.float64, device=self.device)
val[..., 0] = bm.array(self.ux(x, y), device=self.device)
val[..., 1] = bm.array(self.uy(x, y), device=self.device)
return val

@cartesian
Expand All @@ -123,10 +126,12 @@ def hessian(self, p):
#sin = bm.sin
#pi = bm.pi
#cos = bm.cos
val = bm.zeros(p.shape[:-1]+(3, ), dtype=bm.float64)
val[..., 0] = self.uxx(x, y)
val[..., 1] = self.uxy(x, y)
val[..., 2] = self.uyy(x, y)
x = bm.to_numpy(x)
y = bm.to_numpy(y)
val = bm.zeros(p.shape[:-1]+(3, ), dtype=bm.float64, device=self.device)
val[..., 0] = bm.array(self.uxx(x, y), device=self.device)
val[..., 1] = bm.array(self.uxy(x, y), device=self.device)
val[..., 2] = bm.array(self.uyy(x, y), device=self.device)
return val

def dirichlet(self, p):
Expand Down Expand Up @@ -236,7 +241,7 @@ def dirichlet(self, p):
return self.solution(p)


def get_flist(u_sp):
def get_flist(u_sp, device=None):
x = sp.symbols("x")
y = sp.symbols("y")

Expand Down Expand Up @@ -273,40 +278,54 @@ def get_flist(u_sp):
uyyyx = sp.lambdify(('x', 'y'), uyyyx_sp, 'numpy')
uyyyy = sp.lambdify(('x', 'y'), uyyyy_sp, 'numpy')

f = lambda node : u(node[..., 0], node[..., 1])
#f = lambda node : u(node[..., 0], node[..., 1])
def f(node):
x = node[..., 0]
y = node[..., 1]
x_cpu = bm.to_numpy(x)
y_cpu = bm.to_numpy(y)
return bm.array(u(x_cpu, y_cpu), device=device)
def grad_f(node):
x = node[..., 0]
y = node[..., 1]
val = bm.zeros_like(node)
val[..., 0] = ux(x, y)
val[..., 1] = uy(x, y)
val = bm.zeros_like(node, device=device)
x_cpu = bm.to_numpy(x)
y_cpu = bm.to_numpy(y)
val[..., 0] = bm.array(ux(x_cpu, y_cpu), device=device)
val[..., 1] = bm.array(uy(x_cpu, y_cpu), device=device)
return val
def grad_2_f(node):
x = node[..., 0]
y = node[..., 1]
val = bm.zeros(x.shape+(3, ), dtype=bm.float64)
val[..., 0] = uxx(x, y)
val[..., 1] = uyx(x, y)
val[..., 2] = uyy(x, y)
val = bm.zeros(x.shape+(3, ), dtype=bm.float64, device=device)
x = bm.to_numpy(x)
y = bm.to_numpy(y)
val[..., 0] = bm.array(uxx(x, y), device=device)
val[..., 1] = bm.array(uyx(x, y), device=device)
val[..., 2] = bm.array(uyy(x, y), device=device)
return val
def grad_3_f(node):
x = node[..., 0]
y = node[..., 1]
val = bm.zeros(x.shape+(4, ), dtype=bm.float64)
val[..., 0] = uxxx(x, y)
val[..., 1] = uyxx(x, y)
val[..., 2] = uyyx(x, y)
val[..., 3] = uyyy(x, y)
x = bm.to_numpy(x)
y = bm.to_numpy(y)
val = bm.zeros(x.shape+(4, ), dtype=bm.float64, device=device)
val[..., 0] = bm.array(uxxx(x, y), device=device)
val[..., 1] = bm.array(uyxx(x, y), device=device)
val[..., 2] = bm.array(uyyx(x, y), device=device)
val[..., 3] = bm.array(uyyy(x, y), device=device)
return val
def grad_4_f(node):
x = node[..., 0]
y = node[..., 1]
val = bm.zeros(x.shape+(5, ), dtype=bm.float64)
val[..., 0] = uxxxx(x, y)
val[..., 1] = uyxxx(x, y)
val[..., 2] = uyyxx(x, y)
val[..., 3] = uyyyx(x, y)
val[..., 4] = uyyyy(x, y)
x = bm.to_numpy(x)
y = bm.to_numpy(y)
val = bm.zeros(x.shape+(5, ), dtype=bm.float64, device=device)
val[..., 0] = bm.array(uxxxx(x, y), device=device)
val[..., 1] = bm.array(uyxxx(x, y), device=device)
val[..., 2] = bm.array(uyyxx(x, y), device=device)
val[..., 3] = bm.array(uyyyx(x, y), device=device)
val[..., 4] = bm.array(uyyyy(x, y), device=device)
return val

flist = [f, grad_f, grad_2_f, grad_3_f, grad_4_f]
Expand Down
Loading