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

Question: material function with an unexpected keyword 'mode' #623

Closed
Gavinlee960312 opened this issue Jun 23, 2020 · 5 comments
Closed
Labels

Comments

@Gavinlee960312
Copy link

I'm trying to use sfepy to build a mode solver for step_index fiber. But when I tried to define the permittivity in the materials, I cannot find any way to define a circle-shape parameter for materials. Here is my code:

from future import absolute_import
from sfepy import data_dir
import numpy as np
filename_mesh = data_dir + '/meshes/2d/test1.mesh'

lamda=1.55E-6
k0=2*np.pi/lamda
neigen=10
regions={
'Omega':'all',
'Outsider':('vertices of group 1','facet'),
}

fields={
'ElectricField':('complex',2,'Omega',1),
}

variables={
'Ef':('unknown field','ElectricField',0),
'et':('test field','ElectricField','Ef'),
}

ebcs={
'freespace':('Outsider',{'Ef.0':0.0}),
}

materials={
'epsl':'permitivity',
'one':({'one':1.0},),
}

equations={
'ModeEquation':
"""dw_laplace.2.Omega(one.one,et,Ef)
=%sdw_volume_dot.2.Omega(epsl.f,et,Ef)
-%s
dw_volume_dot.2.Omega(one.one,et,Ef)"""
%(k02, neigen*k02)
}

options = {
'nls' : 'newton',
'ls' : 'ls',
}

solvers={
'ls' : ('ls.scipy_direct', {}),
'newton':('nls.newton',{
'i_max':1,
'eps_a':1E-10,
}),
}

def permitivity(coors, mode=None):
if mode=='qp':
r = np.sqrt(coors[:,0]**2.0 + coors[:,1]2.0)
val = 1.45
(r<=60) + 1.4
(r>60 and r <200)
val.shape=(coors.shape[0],1,1)
return {'f':val}
functions = {
'permitivity' : (permitivity,),
}

And the log is:

sfepy: left over: ['name', 'doc', 'package', 'loader', 'spec', 'file', 'cached', 'builtins', 'absolute_import', 'data_dir', 'np', 'lamda', 'k0', 'neigen', 'permitivity', 'verbose', '_filename']
sfepy: reading mesh (C:\Users\Gavin\OneDrive\Desktop\sfepy\meshes\2d\test1.mesh)...
sfepy: number of vertices: 147142
sfepy: number of cells:
sfepy: 2_3: 292260
sfepy: ...done in 1.33 s
sfepy: warning: bad element orientation, trying to correct...
sfepy: ...corrected
sfepy: creating regions...
sfepy: Omega
sfepy: Outsider
sfepy: ...done in 0.18 s
sfepy: equation "ModeEquation":
sfepy: dw_laplace.2.Omega(one.one,et,Ef)
=16432223768723.176dw_volume_dot.2.Omega(epsl.f,et,Ef)
-164322237687231.75
dw_volume_dot.2.Omega(one.one,et,Ef)
sfepy: using solvers:
ts: no ts
nls: newton
ls: ls
sfepy: updating variables...
sfepy: ...done
sfepy: setting up dof connectivities...
sfepy: ...done in 0.03 s
sfepy: matrix shape: (148020, 148020)
sfepy: assembling matrix graph...
sfepy: ...done in 0.10 s
sfepy: matrix structural nonzeros: 1041232 (4.75e-05% fill)
sfepy: updating variables...
sfepy: ...done
sfepy: updating materials...
sfepy: one
sfepy: epsl
Traceback (most recent call last):
File "./simple.py", line 182, in
main()
File "./simple.py", line 179, in main
app()
File "C:\Users\Gavin\OneDrive\Desktop\sfepy\sfepy\applications\application.py", line 29, in call_basic
return self.call(**kwargs)
File "C:\Users\Gavin\OneDrive\Desktop\sfepy\sfepy\applications\pde_solver_app.py", line 228, in call
post_process_hook_final=self.post_process_hook_final)
File "C:\Users\Gavin\OneDrive\Desktop\sfepy\sfepy\discrete\problem.py", line 1435, in solve
status=status)
File "C:\Users\Gavin\OneDrive\Desktop\sfepy\sfepy\solvers\ts_solvers.py", line 34, in _standard_ts_call
status=status, **kwargs)
File "C:\Users\Gavin\OneDrive\Desktop\sfepy\sfepy\solvers\ts_solvers.py", line 71, in call
prestep_fun(ts, vec0)
File "C:\Users\Gavin\OneDrive\Desktop\sfepy\sfepy\discrete\problem.py", line 1260, in prestep_fun
self.update_materials(verbose=self.conf.get('verbose', True))
File "C:\Users\Gavin\OneDrive\Desktop\sfepy\sfepy\discrete\problem.py", line 591, in update_materials
problem=self, verbose=verbose)
File "C:\Users\Gavin\OneDrive\Desktop\sfepy\sfepy\discrete\equations.py", line 340, in time_update_materials
verbose=verbose)
File "C:\Users\Gavin\OneDrive\Desktop\sfepy\sfepy\discrete\materials.py", line 58, in time_update
mat.time_update(ts, equations, mode=mode, problem=problem)
File "C:\Users\Gavin\OneDrive\Desktop\sfepy\sfepy\discrete\materials.py", line 303, in time_update
self.update_data(key, ts, equations, term, problem=problem)
File "C:\Users\Gavin\OneDrive\Desktop\sfepy\sfepy\discrete\materials.py", line 220, in update_data
**self.extra_args)
File "C:\Users\Gavin\OneDrive\Desktop\sfepy\sfepy\discrete\functions.py", line 37, in call
return self.function(*args, **_kwargs)
TypeError: permitivity() got multiple values for argument 'mode'

Anyone can help me? I'll really appreciate it!

@rc
Copy link
Member

rc commented Jun 23, 2020

The material definition functions need to have the following arguments (see https://sfepy.org/doc-devel/users_guide.html#defining-material-parameters-by-functions):

def get_pars(ts, coors, mode=None, **kwargs):
    if mode == 'qp':
        ...

so try using:

def permitivity(ts, coors, mode=None, **kwargs):
    if mode=='qp':
        r = np.sqrt(coors[:,0]**2.0 + coors[:,1]2.0)
        val = 1.45(r<=60) + 1.4(r>60 and r <200)
        val.shape=(coors.shape[0],1,1)
        return {'f':val}

@Gavinlee960312
Copy link
Author

Thanks. I changed the code as you mentioned, but it pops another error:

ValueError: wrong arguments shapes for "+1.0 * dw_laplace.2.Omega(one.one, et, Ef)" term! (see above)

To avoid using dw_laplace, I used dw_div_grad instead, and I got the following error:

File "C:\Users\Gavin\OneDrive\Desktop\sfepy\sfepy\terms\terms_dot.py", line 50, in dw_dot
status = fun(out, mat, val_qp, vgeo, sgeo, fmode)
File "sfepy\terms\extmods\terms.pyx", line 1528, in sfepy.terms.extmods.terms.dw_volume_dot_vector
array2fmfield4(_coef, coef)
File "sfepy\discrete\common\extmods_fmfield.pyx", line 5, in sfepy.discrete.common.extmods._fmfield.array2fmfield4
cdef inline int array2fmfield4(FMField *out,
ValueError: Buffer has wrong number of dimensions (expected 4, got 3)

@rc
Copy link
Member

rc commented Jun 23, 2020

You can use just dw_laplace.2.Omega(et,Ef), the same with dw_volume_dot. Or do you want to use other values in future?

@rc
Copy link
Member

rc commented Jun 23, 2020

Ah, you variable is a vector. So yes, you should use dw_div_grad instead of dw_laplace.

@rc rc added the question label Nov 6, 2020
@rc
Copy link
Member

rc commented Nov 6, 2020

The question was answered, closing.

@rc rc closed this as completed Nov 6, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

2 participants