MPh-py / MPh

Pythonic scripting interface for Comsol Multiphysics
https://mph.readthedocs.io
MIT License
265 stars 67 forks source link

Not all time steps returned from parametric sweep #112

Closed wacht02 closed 1 year ago

wacht02 commented 1 year ago

Greetings. I think i have run into the same problem as #104 where the model.evaluate() function does not return the full solution.

I want to perform a parametric sweep where the solver will take a different amount of time-steps for each parameter. To reproduce this behaviour I modified the capacitor.mph file so that it reflects what I want by making the time-step in the output times range dependent on the parameter d. I also adjusted the Steps taken by the solver selection to strict. This way the solver is forced to take different time-steps for each parameter. The modified file can be found here https://uni-muenster.sciebo.de/s/z8DxEdVFnBRINVF. If you look into the parametric solutions after solving the model you will find, that the parameter d = 1mm has 106 time-steps, the parameter d = 2mm has 205 time-steps and the parameter d = 3mm has 304 time-steps stored in the solution respectivly. When I now open the solved model in python and try to extract the time-steps for the different parameters like with this code:

model = client.load('capacitor_solved.mph') (indices, values) = model.outer('parametric sweep') for i in indices: t = model.evaluate('t', 's', 'parametric sweep', outer = i) print(f'Parameter Value: {values[i-1]} \t Number of Time-steps: {len(t)} \t Last Time: {t[-1]}') client.remove(model)

I get the following print:

Parameter Value: 1.0 Number of Time-steps: 106 Last Time: 1.0 Parameter Value: 2.0 Number of Time-steps: 106 Last Time: 0.505 Parameter Value: 3.0 Number of Time-steps: 106 Last Time: 0.34

So even though the solution for the parameter has more times stored model.evaluate('t', 's', 'parametric sweep') will only return the first 106 time-steps. I think model.evaluate() is somehow limited to the amount of time-steps that are stored in the first parameter value. I checked this by changing the order of the parameters in the sweep. Note that you have to change the order in COMSOL and not in python by i.e. reversing the indices list.

The easy work around for this would be to start with the parameter with the largest number of time-steps stored. But at least for my simulations it is not trivial which parameter that would be and so this option is limited. Another option would be to store the parametric solutions seperatly from each other. Also not realy useable for my problem.

I hope this explains the problem I am having to a sufficient degree. If you already know a solution to this please let me know.

john-hen commented 1 year ago

Hi, thanks a lot for reporting this. I don't know a solution yet, would have to look into this. Not sure when I'll find time for that. But I can reproduce the error, and that is very helpful.

To keep all in one place, I've downloaded the demo model you provided, removed the solutions, and attached it here as time-steps-missing.zip (It's actually the .mph file, so needs to be renamed. GitHub won't accept .mph files as attachments.)

This is the reproducer:

import mph

client = mph.start()
model = client.load('time-steps-missing.mph')
model.solve('sweep')
(indices, values) = model.outer('parametric sweep')
for (index, value) in zip(indices, values):
    times = model.evaluate('t', 's', 'parametric sweep', outer=index)
    print(f'd = {value} mm: {len(times)} time steps')

This the output:

$ python time-steps-missing.py
d = 1.0 mm: 106 time steps
d = 2.0 mm: 106 time steps
d = 3.0 mm: 106 time steps
john-hen commented 1 year ago

Fixed in MPh 1.2.2, released today.

wacht02 commented 1 year ago

Hi, thanks for the quick solution to this problem. I have however run into some more problems that might be linked to this.

While the solution you provided does work for the time expression it doesn't work for some others. To show this behaviour I use the same modified capacitor model as above. With the solved model and using this script

model = client.load('capacitor_solved.mph')
i = 0
exp = ['t', 'ec.Ex']
unit = ['ns', 'V/m']
for (index, value) in zip(indices, values):
    times = model.evaluate(exp[i], unit[i], 'parametric sweep', outer=index)
    print(f'd = {value} mm: {len(times)} time steps')
ende = time.time()
client.remove(model)

to print out the number of time steps for different expressions, I recieve the following result. With i=0 we look at the expression for the time in the model:

d = 1.0 mm: 106 time steps
d = 2.0 mm: 205 time steps
d = 3.0 mm: 304 time steps

And with i = 1 we evaluate the x-compontent of the electric field for all times:

d = 1.0 mm: 106 time steps
d = 2.0 mm: 106 time steps
d = 3.0 mm: 106 time steps

As you can see the first evaluation returns the expected number of timesteps, while the second evaluation only returns a maximum number of timesteps equal to the first outer solution evaluated. So the number of time steps returned depends on the expression i want to evaluate. I am not sure what is causing this issue so if you know a solution please let me know.

john-hen commented 1 year ago

Yeah, okay, that's also wrong.

It is technically a different issue as these are different types of "evaluations". The time t is a global parameter, so it's a "global evaluation" in Comsol terminology. Whereas the electric field is, well, a field, so what Comsol calls a "local evaluation". The returned array then has an extra dimension, for the location. (The coordinates x, y, and z, if this were 3d, would evaluate to one-dimensional arrays that correspond to the same points in space, per index, as the field.)

The previous fix only affects the global expressions. As I've remarked in the commit message, I don't really get why the old way of doing this did not work for your use case. It could be a bug in Comsol. Or the API documentation isn't very clear on what's the right Java method to call. I don't know. But I basically just did it a different way, and that gave the expected result. Though as I understand it, both ways should return the same result, in theory.

I'll have to see if there's a similar work-around for local expressions (fields). I hope there is. But I don't know at this point. Again, I don't see a good demonstration of the "canonical" way to do this kind of thing in the Comsol documentation (for Java), but I may be missing something.

john-hen commented 1 year ago

For what it's worth, I've always avoided setting up parameter sweeps in Comsol if I was automating the model anyway. As then you can just sweep the parameter(s) in Python. Means you'd have to do the entire post-processing in Python though (plots, etc.) as you won't get "one Comsol model with all the solutions" in the end. But my goal was always to discard the solution data, as these are just intermediate results.