Integral over product space

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: 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 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

Okay, I found the problem.
It seems that the potency operator is not thread safe! Investigating the coefficient function of 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 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 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