Dear NGSolve community,
I am having the (maybe little unusual) integral expression in a weak formulation of a problem I am interested in:
[tex]
\int_\Omega \int_\Omega f(x,y,\phi(y)) u(x) v(x) d^3y d^3x.
[/tex]
It is part of a system of equations of a time-dependent problem and both \phi and u are trial functions and v is a test function. I employ a semi-implicit Euler-scheme with Newton’s methods for time integration.
To make things a little easier, I decided to take \phi explicitely, but still the only method I can think of
is to encode an integration scheme for the inner integral in the symbolic expression I give the symbolic integrator in NGSolve. I can imagine that such odd situations are not accounted for right now in NGSolve,
but still wanted to make sure I am not missing a short workaround.
I appreciate any hint, what I might do instead.
All the best,
Philipp
If you take phi explicitly you can just use the gridfunction instead of the trialfunction and integrate it to get the coefficient in your form or am I missing something?
Best
Christopher
Thanks for your fast reply. Since f depends on both x and y, I need the inner integral to be a coefficient function in x. I don’t know a way to do this, but symbolically write down the quadrature in y.
Maybe there is confusion because I named the points x and y which collides with the
notation in NGSolve, where x and y are coordinates. I am sorry for this. Maybe I should edit the post.
All the best,
Philipp
Ah now I got you, sry.
I guess this is an example where you would likely need to write a c++ coefficientfunction to have it somehow performant.
You can use the tutorial here:
https://ngsolve.org/docu/latest/mylittlengs/1_Basic/1_myFEM/coefficient.html
from mylittlengsolve (My CoefficientFunction — NGS-Py 6.2.2302 documentation)
In the evaluate you need to integrate f over the mesh with the fixed x from the mip that you get in the evaluate.
Hope this helps, if you have questions or run into problems with the cpp extension let us know.
Best
I got to work on your suggestion, but sadly got stuck on the following error
that occurs when I import myngspy that I build from the code of “mylittlengsolve”:
ImportError: generic_type: type "MyCoefficient" referenced unknown base type "ngfem::CoefficientFunction"
I modified the example to only implement the MyCoefficient class, so my myngspy.cpp file is slightly different (attached), but that shouldn’t be the problem as far as I understand.
After some digging I suspect that it is related to the pybind-export code for CoefficientFunction not
being executed before my import, but I can’t figure out how I enforce this.
Do you have any suggestions?
Thank you very much
Philipp
Attachment: myngspy.cpp
Hi,
usually the statement py::module::import(“ngsolve”); sould execute the ngsolve python export.
Another option would be to do the “import ngsolve” before importing your module.
Can you zip your whole example and attach it?
Hi,
that is exactly what I was expecting, thus the confusion.
The zip is attached. I kept all paths as in my configuration and
commented out a “fix” (namely explicitely binding CoefficientFunction)
that shouldn’t be necessary.
To my further confusion, importing ngsolve in python doesn’t help.
All the best,
Philipp
Attachment: test_2020-05-26.zip
Hi, after commenting out your stuff in the CMakeLists.txt it worked for me.
Do you use a self compiled ngsolve or the installer?
Do you mean the CoefficientFunction export? With this, it works for me too, but this should be done with the ngsolve import already, right? It was just a “dirty” fix of mine and was not part of the “mylittlengsolve” example.
I compiled my own version.
It is okay for me to export the ngsolve core classes manually. I just thought this might later cause problems since I got something wrong.
No, what I mean I just commented out the setting of the python executable with your hard paths and so (this should actually automatically be set correctly with the find_package(NGSolve))
You should only need this in your CMakeLists.txt:
project(mylittlengsolve)
cmake_minimum_required(VERSION 3.1)
find_package(NGSolve CONFIG REQUIRED
HINTS path_to_ngsolve_dir)
add_ngsolve_python_module(myngspy myngspy.cpp
1_myFEM/myCoefficient.cpp
)
install(TARGETS myngspy DESTINATION .)
There should be no need to export any of the ngsolve core functionality
Ah, sorry. It’s been a long day. You already stated that quite clearly.
I have to check this. I added the python paths
since my system is on python 3.8, but my ngsolve is based on python 3.7.
The python paths I specify are all correct (I just checked). Removing them causes me another error
terminate called after throwing an instance of 'pybind11::error_already_set'
what(): SystemError: <built-in method __contains__ of dict object at 0x7fccec96fc80> returned a result with an error set
Caught SIGABRT: usually caused by abort() or assert()
I think this is highly related to my setup and so I don’t want to bother you with this anymore. Can you maybe just tell me what python path(s) would be responsible for the correct loading of the pybind11 code for the ngsolve internal classes?
Thank you very much for your time!
Philipp
If you use python 3.8 you should build ngsolve against 3.8 as well, since ngsolve/pybind needs to link against the correct python libraries.
I guess you have 2 python versions around and ngsolve and your module find different ones which messes everything up 
Thanks again. You were right. I compiled ngsolve against my globally installed python3.8 and now the error is gone. Sadly I couldn’t find out what exactly was missing for the locally installed python3.7 version, but sometimes it seems better to leave it that way ;).
All the best,
Philipp
I managed to implement a coefficient function that represents
[tex]
(xx, yy) \mapsto \int_\Omega f((xx,yy),(x,y)) d(x,y);
[/tex] the code is attached. Single threaded it works quite well, but I get odd NaNs at some points
(xx,yy) when computing
[tex]
\int_{[0,1]^2} \sqrt{ (xx - 0.5)^2 + (yy - 0.5)^2 } d(x,y)
[/tex]
with the TaskManager. The following python code should reproduce the error:
from ngsolve import *
from netgen.geom2d import unit_square
from intcoefffun import XArgumentMarker, YArgumentMarker, IntCoefficientFunction
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
xx = XArgumentMarker()
yy = YArgumentMarker()
fun = lambda c1, c2: sqrt((c1-0.5)**2 + (c2-0.5)**2)
intcoeff = IntCoefficientFunction(mesh, fun(xx,yy))
with TaskManager():
print("int = ", Integrate(intcoeff, mesh))
My question is: Did I made some obvious mistake in my implementation?
So far, I have not been able to narrow down the problem more.
Edit:
I did some more debugging and found that changing
fun = lambda c1, c2: CoefficientFunction(c1-0.5),(c2-0.5)).Norm()
works well. Glad that I could work around this problem I went on just
to discover that I now get other race conditions with coefficient
functions as exp and acos. They show in integral values varying
in a large range. Any idea how this might be related?
All the best,
Philipp
Attachment: intcoeff.zip
Okay, I found the problem.
It seems that the potency operator
x**2
is not thread safe! Investigating the coefficient function of
sqrt((xx-0.5)**2 + (yy-0.5)**2)
in gdb I found that the nodes in which (xx-0.5)**2 and (yy-0.5)**2 should be placed are
not initialised. However, employing pow like
sqrt(pow(xx-0.5, 2) + pow(yy-0.5, 2))
works!
Maybe this is a problem specific to my system, but I want to mention it anyhow…maybe a desperate soul runs into this thread with a similar problem.
Regards,
Philipp
As life is not so simple, I experience some other data race issues with the special coefficient function IfPos.
After some general considerations, I came up with the following question:
Having two integration procedures both acessing the ParameterCoefficientFunctions
XArgumentMarker and YArgumentMarker, do I have to introduce a lock?
I mean the “outer” integration procedure triggered by
Integrate(intcoeff, mesh)
alters XArgumentMarker and YArgumentMarker used by the integration
in the coefficient function possibly not waiting until the inner integration is finished.
(I am not familiar with the evaluations going on in the ngsolve core and if they are synchronised.)
If so, I would just use the pthread locking mechanisms, but is there a recommended NGSolve solution for this?
All the best,
Philipp
Hi Philipp,
you found the problem !
Yes, setting the outer coordinates when calling the IntCoefficientFunction to the Parameters changes the expression tree, and does now allow reentrant evaluation.
Of course, you can lock … But then you have to ask what scaling you get, and why to use parallel evaluation (for the outer loop) anyway.
We had similar issues to get element-specific data down to the CoefficientFunction:
Where to store that data ? In specific cases (it should work out for you) one can use thread-local variables to store the outer x,y coordinates.
We have attached the thread-local data to the only variable passed into the Evaluate-Function, the MappedIntegrationPoint/Rule. It points to the ElementTransformation, to which we can assign user-specific data to the void-pointer userdata. See use-cases in fem/symbolicintegrator.cpp
You implement Evaluate - methods in the ArgumentMarker classes, which unpack the userdata.
I recommend that approach.
Anyway, it looks like you have already dug quite well into NGSolve-C++
Best,
Joachim
Hi Joachim,
thank you very much for the assurance and the detailed explanation. That sounds pretty good. I
thought about copying the coefficient function’s expression tree before iterating
over the mesh in the inner integration, but your solutions seems much better. It stores less data
and since the structure for thread-specific data is already there, I shall use it.
I will try tomorrow.
All the best,
Philipp
After some testing, I’ve convinced myself of my implementation being sensible.
I attached the code in case anybody else might be interested.
Thank to Joachim and Christopher, again!
All the best,
Philipp
Attachment: integral_coefficient_function.zip