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.

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

Non Matching Multipatch #320

Merged
merged 81 commits into from
Aug 2, 2024
Merged
Changes from 1 commit
Commits
Show all changes
81 commits
Select commit Hold shift + click to select a range
d830ae1
add non matching conforming projection operators
FrederikSchnack Jul 14, 2023
a27a062
add some plotting utilities
FrederikSchnack Jul 14, 2023
bbad569
make multipatch examples run
FrederikSchnack Jul 14, 2023
ec10dfd
Use SymPDE development version from GitHub
yguclu Jul 17, 2023
a166d02
allow for non-matching grids in postprocessing
FrederikSchnack Jul 17, 2023
87bbef6
add utils
FrederikSchnack Jul 17, 2023
1a41cd0
adds non-matching domain utilities
FrederikSchnack Jul 18, 2023
e84552a
adds curl-curl test case comparisson
FrederikSchnack Jul 18, 2023
548881a
small changes
FrederikSchnack Jul 24, 2023
b0e64c9
Merge remote-tracking branch 'origin/devel' into non_matching_multipatch
FrederikSchnack Jul 24, 2023
4dbc422
Merge branch 'devel' into non_matching_multipatch
yguclu Jul 25, 2023
b3d1451
add time-domain Maxwell
FrederikSchnack Jul 27, 2023
a806b6d
add time harmonic and time domain maxwell examples
FrederikSchnack Jul 28, 2023
26578f6
add new tests and last example
FrederikSchnack Aug 2, 2023
f68c562
Merge branch 'devel' into non_matching_multipatch
FrederikSchnack Sep 20, 2023
7fbaecf
make Codacy happy
FrederikSchnack Sep 20, 2023
b52041f
add pml example
FrederikSchnack Nov 27, 2023
7f5be32
update non_matching operators and add proposal of tests
FrederikSchnack Nov 27, 2023
1692570
adapt files to new projection
FrederikSchnack Nov 27, 2023
544a237
make codacy happy
FrederikSchnack Nov 28, 2023
7b357bb
add some pml experiments
FrederikSchnack Dec 12, 2023
ebaa9f3
Merge branch 'devel' into non_matching_multipatch
FrederikSchnack Mar 1, 2024
a147990
add docs and readability
FrederikSchnack Mar 5, 2024
c7325da
Merge branch 'devel' into non_matching_multipatch
FrederikSchnack Mar 5, 2024
567ce21
move pml experiments to its own branch
FrederikSchnack Mar 5, 2024
2af29a8
adapt tests
FrederikSchnack Mar 5, 2024
d4fbdbf
change Hodge matrix naming conventions
FrederikSchnack Mar 7, 2024
11a9859
Merge branch 'devel' into non_matching_multipatch
FrederikSchnack Mar 7, 2024
ae90d77
Merge branch 'devel' into non_matching_multipatch
yguclu Mar 19, 2024
05bd083
Merge branch 'devel' into non_matching_multipatch
FrederikSchnack Apr 23, 2024
ff8c1c0
adds revised Hcurl conforming projections and adapts the tests, prior…
FrederikSchnack May 22, 2024
53c1b3c
adapt conforming examples to new projections
FrederikSchnack May 22, 2024
dcc86b1
update non-matching examples
FrederikSchnack May 23, 2024
a853a47
work in comments
FrederikSchnack May 23, 2024
5371468
update tests
FrederikSchnack May 23, 2024
9357bdd
Merge branch 'devel' into non_matching_multipatch
FrederikSchnack May 23, 2024
ba230ba
stop using create_domain() -- wip
campospinto May 31, 2024
685f4a5
use Domain.join() to build non-matching domains
campospinto Jun 2, 2024
a3a5341
using temp function for sympde Domain.join
campospinto Jun 4, 2024
220f1cf
adding ref
campospinto Jun 4, 2024
52fec60
Merge example scripts and get rid of get_source_and_solution_OBSOLETE
FrederikSchnack Jun 10, 2024
4d5f780
clean up multipatch domain utilities
FrederikSchnack Jun 11, 2024
fffe43f
Coarsen test runs and add timedomain dummy run
FrederikSchnack Jun 11, 2024
d417da4
Make codacy happy
FrederikSchnack Jun 12, 2024
a748a4d
fix typo
FrederikSchnack Jun 12, 2024
2bceaf6
Expose multipatch modules to docs
FrederikSchnack Jun 17, 2024
383d341
minor documentation related fixes
jowezarek Jun 17, 2024
a9da50e
sphinx fix
campospinto Jun 18, 2024
7f9b522
functions stencil index to petsc index and viceversa
e-moral-sanchez May 8, 2024
b8bdab9
efficient conversion to PETSc of StencilVectors
e-moral-sanchez May 8, 2024
eba63ff
works for 1D stencil matrix with multiple processes
e-moral-sanchez May 13, 2024
20b7fe7
correct indexing for 2D stencilmatrix
e-moral-sanchez May 15, 2024
f11599d
works 1d,2d,3d stencilmatrix without periodic BC
e-moral-sanchez May 15, 2024
7a2c555
efficient conversion of stencilmatrix to PETSc.Mat for 1,2,3D and pe…
e-moral-sanchez May 16, 2024
4ad477c
Fixed conversion from 2D BlockStencilVector to Petsc.Vec
e-moral-sanchez May 17, 2024
f62385e
fix general case also for stencilvector 2D
e-moral-sanchez May 17, 2024
1256836
fixed petsc_to_psydac for BlockVectors
e-moral-sanchez May 21, 2024
b228eac
PETSc conversion works for StencilVector and BlockStencilVector of 1 …
e-moral-sanchez May 21, 2024
437c38f
conversion works for StencilMatrix
e-moral-sanchez May 21, 2024
80f3e6c
works for BlockLinearOperators, the blocks of which are Stencilmatrices
e-moral-sanchez May 21, 2024
abe2614
Clean up, docstrings
e-moral-sanchez May 21, 2024
21f06dc
fix bugs
e-moral-sanchez May 22, 2024
889619d
sequential case, docstrings
e-moral-sanchez May 22, 2024
ed5abc0
cleaning
e-moral-sanchez May 22, 2024
cddbd38
erase forgotten comments
e-moral-sanchez May 22, 2024
bd2fcd1
fix loop stencil
e-moral-sanchez May 23, 2024
ab228b9
forgotten prints
e-moral-sanchez May 27, 2024
820bd7b
add pyccel kernel
e-moral-sanchez May 29, 2024
3a7c2a7
string for type annotations
e-moral-sanchez May 29, 2024
1778055
clean up plus fix serial case
e-moral-sanchez May 29, 2024
bdc0645
clean time measurements
e-moral-sanchez Jun 3, 2024
a560b90
fix version 12 of macos in tests
e-moral-sanchez Jun 12, 2024
f65b21e
Fixes for CI failures caused by new macOS runner version and numpy 2.0
kvrigor Jun 19, 2024
9172dfc
CI: Temporarily disabled docs deployment since base repo is a fork
kvrigor Jun 19, 2024
45cf90c
Update macos version
e-moral-sanchez Jun 21, 2024
ce4f189
Disable deploy_docs in documentation.yml
FrederikSchnack Jun 21, 2024
6ff11ee
Revert "Disable deploy_docs in documentation.yml"
yguclu Aug 1, 2024
613b077
Revert "CI: Temporarily disabled docs deployment since base repo is a…
yguclu Aug 1, 2024
bc75682
Style cleanup of stencil2IJV_kernels.py
yguclu Jul 31, 2024
8d429c1
Apply some PEP8 recommendations to topetsc.py
yguclu Aug 1, 2024
eaa9854
Merge branch 'devel' into non_matching_multipatch
yguclu Aug 1, 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
Prev Previous commit
Next Next commit
add time-domain Maxwell
  • Loading branch information
FrederikSchnack committed Jul 27, 2023
commit b3d1451cb19a93ad1ffc4ec345a27238e9bdc119
340 changes: 292 additions & 48 deletions psydac/feec/multipatch/examples/ppc_test_cases.py
Original file line number Diff line number Diff line change
@@ -5,7 +5,7 @@
import os
import numpy as np

from sympy import pi, cos, sin, Tuple, exp
from sympy import pi, cos, sin, Tuple, exp, atan, atan2

from sympde.topology import Derham

@@ -16,18 +16,145 @@
from psydac.feec.multipatch.plotting_utilities import get_plotting_grid, my_small_plot, my_small_streamplot
from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain

from psydac.feec.multipatch.utilities import sol_ref_fn, error_fn, get_method_name, get_fem_name, get_load_dir

comm = MPI.COMM_WORLD


# todo [MCP, 12/02/2022]: add an 'equation' argument to be able to return 'exact solution'

def get_phi_pulse(x_0, y_0, domain=None):
x,y = domain.coordinates
ds2_0 = (0.02)**2
sigma_0 = (x-x_0)**2 + (y-y_0)**2
phi_0 = exp(-sigma_0**2/(2*ds2_0))

return phi_0

def get_div_free_pulse(x_0, y_0, domain=None):
x,y = domain.coordinates
ds2_0 = (0.02)**2
sigma_0 = (x-x_0)**2 + (y-y_0)**2
phi_0 = exp(-sigma_0**2/(2*ds2_0))
dx_sig_0 = 2*(x-x_0)
dy_sig_0 = 2*(y-y_0)
dx_phi_0 = - dx_sig_0 * sigma_0 / ds2_0 * phi_0
dy_phi_0 = - dy_sig_0 * sigma_0 / ds2_0 * phi_0
f_x = dy_phi_0
f_y = - dx_phi_0
f_vect = Tuple(f_x, f_y)

return f_vect

def get_curl_free_pulse(x_0, y_0, domain=None, pp=False):
# return -grad phi_0
x,y = domain.coordinates
if pp:
# psi=phi
ds2_0 = (0.02)**2
else:
ds2_0 = (0.1)**2
sigma_0 = (x-x_0)**2 + (y-y_0)**2
phi_0 = exp(-sigma_0**2/(2*ds2_0))
dx_sig_0 = 2*(x-x_0)
dy_sig_0 = 2*(y-y_0)
dx_phi_0 = - dx_sig_0 * sigma_0 / ds2_0 * phi_0
dy_phi_0 = - dy_sig_0 * sigma_0 / ds2_0 * phi_0
f_x = -dx_phi_0
f_y = -dy_phi_0
f_vect = Tuple(f_x, f_y)

return f_vect

def get_Delta_phi_pulse(x_0, y_0, domain=None, pp=False):
# return -Delta phi_0, with same phi_0 as in get_curl_free_pulse()
x,y = domain.coordinates
if pp:
# psi=phi
ds2_0 = (0.02)**2
else:
ds2_0 = (0.1)**2
sigma_0 = (x-x_0)**2 + (y-y_0)**2
phi_0 = exp(-sigma_0**2/(2*ds2_0))
dx_sig_0 = 2*(x-x_0)
dy_sig_0 = 2*(y-y_0)
dxx_sig_0 = 2
dyy_sig_0 = 2
dxx_phi_0 = ((dx_sig_0 * sigma_0 / ds2_0)**2 - ((dx_sig_0)**2 + dxx_sig_0 * sigma_0)/ds2_0 ) * phi_0
dyy_phi_0 = ((dy_sig_0 * sigma_0 / ds2_0)**2 - ((dy_sig_0)**2 + dyy_sig_0 * sigma_0)/ds2_0 ) * phi_0
f = - dxx_phi_0 - dyy_phi_0

return f

def get_Gaussian_beam(x_0, y_0, domain=None):
# return E = cos(k*x) exp( - x^2 + y^2 / 2 sigma^2) v
x,y = domain.coordinates
x = x - x_0
y = y - y_0

k = 2*pi
sigma = 0.7

xy = x**2 + y**2
ef = exp( - xy/(2*sigma**2) )

E = cos(k * y) * ef
B = y/(sigma**2) * E - sin(k * y) * ef

return Tuple(E, 0), B

def get_Gaussian_beam2(x_0, y_0, domain=None):
# return E = cos(k*x) exp( - x^2 + y^2 / 2 sigma^2) v
x,y = domain.coordinates


x0 = x_0
y0 = y_0
theta = pi/2
w0 = 1

t = [(x-x0)*cos(theta) - (y - y0) * sin(theta), (x-x0)*sin(theta) + (y-y0) * cos(theta)]

## Gaussian beam
'''Beam inciding from the left, centered and normal to wall:
x: axial normalized distance to the beam's focus
y: radial normalized distance to the center axis of the beam
'''
EW0 = 1.0 # amplitude at the waist
k0 = 2 * pi # free-space wavenumber

x_ray = pi * w0 ** 2 # Rayleigh range

w = w0 * ( 1 + t[0]**2/x_ray**2 )**0.5 #width
curv = t[0] / ( t[0]**2 + x_ray**2 ) #curvature

gouy_psi = -0.5 * atan2(t[0] / x_ray, 1.) # corresponds to atan(x / x_ray), which is the Gouy phase

EW_mod = EW0 * (w0 / w)**0.5 * exp(-(t[1] ** 2) / (w ** 2)) # Amplitude
phase = k0 * t[0] + 0.5 * k0 * curv * t[1] ** 2 + gouy_psi # Phase

EW_r = EW_mod * cos(phase) # Real part
EW_i = EW_mod * sin(phase) # Imaginary part

B = 0#t[1]/(w**2) * EW_r

return Tuple(0,EW_r), B


def get_source_and_sol_for_magnetostatic_pbm(
source_type=None,
domain=None, domain_name=None,
refsol_params=None
):
"""
provide source, and exact solutions when available, for:

Find u=B in H(curl) such that

div B = 0
curl B = j

written as a mixed problem, see solve_magnetostatic_pbm()
"""
u_ex = None # exact solution
x,y = domain.coordinates
if source_type == 'dipole_J':
# we compute two possible source terms:
@@ -62,23 +189,171 @@ def get_source_and_sol_for_magnetostatic_pbm(
else:
raise ValueError(source_type)

# ref solution in V1h:
uh_ref = get_sol_ref_V1h(source_type, domain, domain_name, refsol_params)
return f_scal, f_vect, j_scal, u_ex

return f_scal, f_vect, j_scal, uh_ref

def get_source_and_solution_hcurl(
source_type=None, eta=0, mu=0, nu=0,
domain=None, domain_name=None):
"""
provide source, and exact solutions when available, for:

Find u in H(curl) such that

A u = f on \Omega
n x u = n x u_bc on \partial \Omega

with

A u := eta * u + mu * curl curl u - nu * grad div u

see solve_hcurl_source_pbm()
"""

# exact solutions (if available)
u_ex = None
curl_u_ex = None
div_u_ex = None

# bc solution: describe the bc on boundary. Inside domain, values should not matter. Homogeneous bc will be used if None
u_bc = None

# source terms
f_vect = None

# auxiliary term (for more diagnostics)
grad_phi = None
phi = None

x,y = domain.coordinates

if source_type == 'manu_maxwell_inhom':
# used for Maxwell equation with manufactured solution
f_vect = Tuple(eta*sin(pi*y) - pi**2*sin(pi*y)*cos(pi*x) + pi**2*sin(pi*y),
eta*sin(pi*x)*cos(pi*y) + pi**2*sin(pi*x)*cos(pi*y))
if nu == 0:
u_ex = Tuple(sin(pi*y), sin(pi*x)*cos(pi*y))
curl_u_ex = pi*(cos(pi*x)*cos(pi*y) - cos(pi*y))
div_u_ex = -pi*sin(pi*x)*sin(pi*y)
else:
raise NotImplementedError
u_bc = u_ex

elif source_type == 'elliptic_J':
# no manufactured solution for Maxwell pbm
x0 = 1.5
y0 = 1.5
s = (x-x0) - (y-y0)
t = (x-x0) + (y-y0)
a = (1/1.9)**2
b = (1/1.2)**2
sigma2 = 0.0121
tau = a*s**2 + b*t**2 - 1
phi = exp(-tau**2/(2*sigma2))
dx_tau = 2*( a*s + b*t)
dy_tau = 2*(-a*s + b*t)

f_x = dy_tau * phi
f_y = - dx_tau * phi
f_vect = Tuple(f_x, f_y)

else:
raise ValueError(source_type)

def get_source_and_solution(source_type=None, eta=0, mu=0, nu=0,
# u_ex = Tuple(0, 1) # DEBUG
return f_vect, u_bc, u_ex, curl_u_ex, div_u_ex #, phi, grad_phi

def get_source_and_solution_h1(source_type=None, eta=0, mu=0,
domain=None, domain_name=None):
"""
provide source, and exact solutions when available, for:

Find u in H^1, such that

A u = f on \Omega
u = u_bc on \partial \Omega

with

A u := eta * u - mu * div grad u

see solve_h1_source_pbm()
"""

# exact solutions (if available)
u_ex = None

# bc solution: describe the bc on boundary. Inside domain, values should not matter. Homogeneous bc will be used if None
u_bc = None

# source terms
f_scal = None

# auxiliary term (for more diagnostics)
grad_phi = None
phi = None

x,y = domain.coordinates

if source_type in ['manu_poisson_elliptic']:
x0 = 1.5
y0 = 1.5
s = (x-x0) - (y-y0)
t = (x-x0) + (y-y0)
a = (1/1.9)**2
b = (1/1.2)**2
sigma2 = 0.0121
tau = a*s**2 + b*t**2 - 1
phi = exp(-tau**2/(2*sigma2))
dx_tau = 2*( a*s + b*t)
dy_tau = 2*(-a*s + b*t)
dxx_tau = 2*(a + b)
dyy_tau = 2*(a + b)

dx_phi = (-tau*dx_tau/sigma2)*phi
dy_phi = (-tau*dy_tau/sigma2)*phi
grad_phi = Tuple(dx_phi, dy_phi)

f_scal = -( (tau*dx_tau/sigma2)**2 - (tau*dxx_tau + dx_tau**2)/sigma2
+(tau*dy_tau/sigma2)**2 - (tau*dyy_tau + dy_tau**2)/sigma2 )*phi

# exact solution of -p'' = f with hom. bc's on pretzel domain
if mu == 1 and eta == 0:
u_ex = phi
else:
print('WARNING (54375385643): exact solution not available in this case!')

if not domain_name in ['pretzel', 'pretzel_f']:
# we may have non-hom bc's
u_bc = u_ex

elif source_type == 'manu_poisson_2':
f_scal = -4
if mu == 1 and eta == 0:
u_ex = x**2+y**2
else:
raise NotImplementedError
u_bc = u_ex

elif source_type == 'manu_poisson_sincos':
u_ex = sin(pi*x)*cos(pi*y)
f_scal = (eta + 2*mu*pi**2) * u_ex
u_bc = u_ex

else:
raise ValueError(source_type)

return f_scal, u_bc, u_ex



def get_source_and_solution_OBSOLETE(source_type=None, eta=0, mu=0, nu=0,
domain=None, domain_name=None,
refsol_params=None):
"""
compute source and reference solution (exact, or reference values) when possible, depending on the source_type
OBSOLETE: kept for some test-cases
"""
FrederikSchnack marked this conversation as resolved.
Show resolved Hide resolved

# ref solution (values on diag grid)
ph_ref = None
uh_ref = None

# exact solutions (if available)
u_ex = None
p_ex = None
@@ -119,7 +394,7 @@ def get_source_and_solution(source_type=None, eta=0, mu=0, nu=0,

elif source_type == 'manutor_poisson':
# todo: remove if not used ?
# same as manu_poisson, with arbitrary value for tor
# same as manu_poisson_ellip, with arbitrary value for tor
x0 = 1.5
y0 = 1.5
s = (x-x0) - (y-y0)
@@ -174,6 +449,9 @@ def get_source_and_solution(source_type=None, eta=0, mu=0, nu=0,
# exact solution of -p'' = f with hom. bc's on pretzel domain
p_ex = phi

if source_type == 'manu_poisson' and mu == 1 and eta == 0:
u_ex = phi

if not domain_name in ['pretzel', 'pretzel_f']:
print("WARNING (87656547) -- I'm not sure we have an exact solution -- check the bc's on the domain "+domain_name)
# raise NotImplementedError(domain_name)
@@ -325,40 +603,6 @@ def get_source_and_solution(source_type=None, eta=0, mu=0, nu=0,
else:
raise ValueError(source_type)

if u_ex is None:
uh_ref = get_sol_ref_V1h(source_type, domain, domain_name, refsol_params)

return f_scal, f_vect, u_bc, ph_ref, uh_ref, p_ex, u_ex, phi, grad_phi

return f_scal, f_vect, u_bc, p_ex, u_ex, phi, grad_phi

def get_sol_ref_V1h( source_type=None, domain=None, domain_name=None, refsol_params=None ):
"""
get a reference solution as a V1h FemField
"""
uh_ref = None
if refsol_params is not None:
N_diag, method_ref, source_proj_ref = refsol_params
u_ref_filename = ( get_load_dir(method=method_ref, domain_name=domain_name,nc=None,deg=None,data='solutions')
+ sol_ref_fn(source_type, N_diag, source_proj=source_proj_ref) )
print("no exact solution for this test-case, looking for ref solution values in file {}...".format(u_ref_filename))
if os.path.isfile(u_ref_filename):
print("-- file found")
with open(u_ref_filename, 'rb') as file:
ncells_degree = np.load(file)
ncells = [int(i) for i in ncells_degree['ncells_degree'][0]]
degree = [int(i) for i in ncells_degree['ncells_degree'][1]]

derham = Derham(domain, ["H1", "Hcurl", "L2"])
domain_h = discretize(domain, ncells=ncells, comm=comm)
V1h = discretize(derham.V1, domain_h, degree=degree, basis='M')
uh_ref = FemField(V1h)
for i,Vi in enumerate(V1h.spaces):
for j,Vij in enumerate(Vi.spaces):
filename = u_ref_filename+'_%d_%d'%(i,j)
uij = Vij.import_fields(filename, 'phi')
uh_ref.fields[i].fields[j].coeffs._data = uij[0].coeffs._data

else:
print("-- no file, skipping it")

return uh_ref
Loading