Discussion :

[Heimdall]

Salut,

Soit une équation vectorielle S = (1-d^2Laplacien)B.
S et B sont des vecteurs a trois composantes Sx,Sy,Sz et Bx,By, Bz, ne dépendant chacune que de x et y, nous avons donc 3 équations 2D linéaires.

On connait S, d est une constante, et on cherche B.

Les conditions périodiques sont de trois types :

1/ x,y périodiques, on résoud alors les équations dans la base de fourier bien adaptée. Pas de problemes.

2/ x, y non périodiques. Je veux alors des conditions (venant des équations de maxwell) aux bords de Neumann. i.e. je veux que la dérivée normale de la composante normale soit nulle, et que la dérivée normale de la composante tangentielle soit nulle (pour l'instant, on pourra faire ensuite des conditions ou elle n'est pas nulle)

A titre d'exemple, supposons que l'on se place au bord x = L, je veux donc imposer :

a/ dBx/dx = 0
b/ dBy/dx = 0
c/ dBz/dx = 0

Aucun probleme, il suffit d'utiliser la transformée en cosinus, qui assure des conditions de Neumann homogène (dérivée normale = 0 aux bords)

J'ai testé ceci, ça marche bien. Seulement je me suis rendu compte que je n'avais pas div(B) = 0 au bord. Et puis, div(B) = dBx/dx + dBy/dy. Le premier terme est nul par la transformée en cosinus, ok... mais le second (dérivée tangentielle de la composante tangentielle) n'est pas imposé, et donc j'ai div(B) != 0.

Auriez-vous une idée pour assurer également cette condition ?


3/ Conditions limites mixes, périodique en x, et Neuman+div(B)=0 en y. Donc ici, ça correspondrait au cas 1 en x, et cas 2 en y. Seulement là je ne vois pas trop comment faire, idéalement il faudrait faire une TF suivant la direction x et une TC suivant la direction y, mais je me vois mal faire ça avec fftw...


Merci pour vos idées

[FR119492]

Salut !

Les conditions périodiques sont de trois types

des conditions (venant des équations de maxwell) aux bords de Neumann

La notion de périodicité s'applique théoriquement à un domaine s'étendant jusqu'à l'infini, Alors, il n'y a pas de bords.

Je ne suis pas certain, mais je pense que c'est cette histoire de "périodicité" qui embrouille tout.

Jean-Marc Blanc

[Heimdall]

ah je me suis trompé, il faut lire "les conditions limites sont de trois types : ".


En fait j'ai potentiellement une idée. Lorsque je résoud mon équation dans l'espace des cosinus pour less TROIS composantes (au bords : 1 normale et 2 tangentielles) j'impose que la dérivée normale (au bord) des 3 composantes soit nulle.

Plaçons nous au bord x=L par exemple

La résolution en cosinus impose donc :

dBx/dx = dBy/dx = dBz/dx = 0

Et donc pour que mon champ solution soit a divergence nulle, il faut également imposer dBy/dy = 0 (d/dz = 0 est automatiquement satisfait car je suis en 2D)
Ce qui est impossible (je pense) à réaliser avec ma résolution en cosinus.

Une solution serait alors, de faire pour les composantes tangentielle, une résolution en cosinus (condition de Neumann homogène) et une résolution en sinus pour la composante normale (condition de dirichlet homogene).

Cette méthode a l'intéret d'imposer (sur x=L par exemple) :

Bx = 0
dBy/dx = dBz/dx = 0

On ne touche plus aux termes de div(B) dBx/dx et dBy/dy qui satisfont de manière intrinsèque a mon sytsteme d'équation div(B) = 0

Le probleme avant c'est que je touchais (via la condition limite) à un terme de ma divergence sans toucher l'autre, qui était donné par mon code.

Je ne suis pas sûr que ça marche, mais c'est une idée

[FR119492]

Salut !

J'essaie de reformuler ton problème tel qu'il se pose, en élaguant tout ce que tu y as ajouté pour le rendre pratiquement insoluble:

On a un domaine rectangulaire D dans le plan Oxy:
0 =< x =< X
0 =< y =< Y
dans lequel une grandeur vectorielle B, de composantes Bx, By et Bz obéit à l'équation
B - d^2 * laplacien(B) = S
On sait en outre que B est à flux conservatif.
Le coefficient d est une constante scalaire qui a la dimension d'une longueur et S est une grandeur vectorielle donnée, de composantes Sx, Sy et Sz.
On connait aussi les conditions aux limites (indique-nous lesquelles sont justes)
Bx(0,y) = 0 (Dirichlet) ou dBx/dx(0,y) = 0 (Neumann)
By(0,y) = 0 (Dirichlet) ou dBy/dx(0,y) = 0 (Neumann)
By(0,y) = 0 (Dirichlet) ou dBz/dx(0,y) = 0 (Neumann)

Bx(X,y) = 0 (Dirichlet) ou dBx/dx(X,y) = 0 (Neumann)
By(X,y) = 0 (Dirichlet) ou dBy/dx(X,y) = 0 (Neumann)
By(X,y) = 0 (Dirichlet) ou dBz/dx(X,y) = 0 (Neumann)

Bx(x,0) = 0 (Dirichlet) ou dBx/dy(x,0) = 0 (Neumann)
By(x,0) = 0 (Dirichlet) ou dBy/dy(x,0) = 0 (Neumann)
By(x,0) = 0 (Dirichlet) ou dBz/dy(x,0) = 0 (Neumann)

Bx(x,Y) = 0 (Dirichlet) ou dBx/dy(x,Y) = 0 (Neumann)
By(x,Y) = 0 (Dirichlet) ou dBy/dy(x,Y) = 0 (Neumann)
By(x,Y) = 0 (Dirichlet) ou dBz/dy(x,Y) = 0 (Neumann)

Une fois le problème bien posé, la solution sera facile.

Jean-Marc Blanc

[Heimdall]

Salut !

J'essaie de reformuler ton problème tel qu'il se pose, en élaguant tout ce que tu y as ajouté pour le rendre pratiquement insoluble:

On a un domaine rectangulaire D dans le plan Oxy:
0 =< x =< X
0 =< y =< Y
dans lequel une grandeur vectorielle B, de composantes Bx, By et Bz obéit à l'équation
B - d^2 * laplacien(B) = S
On sait en outre que B est à flux conservatif.
Le coefficient d est une constante scalaire qui a la dimension d'une longueur et S est une grandeur vectorielle donnée, de composantes Sx, Sy et Sz.

pour l'instant c'est ça.

On connait aussi les conditions aux limites (indique-nous lesquelles sont justes)

ok je ne laisse que celle que je veux.

bord gauche et droit (x=0 ET x=X, y)
-----------------------------------------------
Bx(0,y) = Bx(X,y) périodique
By(0,y) = By(X,y) périodique
Bz(0,y) = Bz(X,y) périodique


Bord inférieur (x, y=0)
------------------------------
dBx/dy(x,0) = 0 (Neumann)
By(x,0) = 0 (Dirichlet)
dBz/dy(x,0) = 0 (Neumann)


bord supérieur (x, y=Y)
-------------------------------
dBx/dy(x,Y) = 0 (Neumann)
By(x,Y) = 0 (Dirichlet)
dBz/dy(x,Y) = 0 (Neumann)


Si j'ai nx+1 et ny+1 points suivant x et y respectivement, mon idée était de faire la résolution par (ny+1) DFT suivant x, (nx+1) DCT pour (Bx et Bz) suivant y et (nx+1) DST suivant y pour By.

[FR119492]

Salut !

Bx(0,y) = Bx(X,y) périodique
By(0,y) = By(X,y) périodique
Bz(0,y) = Bz(X,y) périodique

Ici, "périodique" ne veut rien dire tant que tu ne précises pas quelle est la période.

La relation peut être satisfaite si la fonction n'est pas périodique.
Bx = x^2 - X*x + X^2/4

La fonction peut être périodique sans que la relation ne soit satisfaite:
Bx = cos(0,5*pi*x/X)

Jean-Marc Blanc

[Heimdall]

Salut !



Ici, "périodique" ne veut rien dire tant que tu ne précises pas quelle est la période.

La relation peut être satisfaite si la fonction n'est pas périodique.
Bx = x^2 - X*x + X^2/4

La fonction peut être périodique sans que la relation ne soit satisfaite:
Bx = cos(0,5*pi*x/X)

Jean-Marc Blanc

Il me semble bien évident que ce que je disais était que la fonction est périodique et de période la taille de la boite. Comment pourrait-on qualifier de périodique une fonction dont la période est plus grande que la boite.

bref.. périodique sur la boite

[Ehouarn]

Bonjour,

Une petite question au sujet de ton système: tu parles de résoudre ton équation (de Helmoltz?) et du fait que B est à divergence nulle. Mais comment couples-tu les deux (il ne suffit pas d'imposer div B=0 seulement aux frontières pour que cela soit vérifié partout dans le domaine, n'est ce pas)?

Ehouarn

[Heimdall]

Bonjour,

Une petite question au sujet de ton système: tu parles de résoudre ton équation (de Helmoltz?) et du fait que B est à divergence nulle. Mais comment couples-tu les deux (il ne suffit pas d'imposer div B=0 seulement aux frontières pour que cela soit vérifié partout dans le domaine, n'est ce pas)?

Ehouarn

Tu as sans doute raison de poser cette question. Imagine un code d'électromagnétisme, ou a chaque pas de temps te résouds E et B, ton champ électrique et ton champ magnétique. Il se trouve que dans mon code, le systeme d'équations implique que div(B)=0 partout. Maintenant j'ai rajouté une pichnette dans le code, qui fait que je résoud non plus E et B, mais E et S, ou S = B - d2*laplacien(B). Il me faut une fois S calculé, retrouver B en inversant ce systeme. S est donné a divergence nulle (c'est le résultat de rot(E)) j'ai vérifié.

J'ai inversé le systeme avec des conditions de Neumann sur tous les bords pour toutes les composantes (qui tente rien n'a rien) et j'ai vu div(B)=0 partout sauf sur les bords. Pourquoi ? Parce que normalement le code assure tout seul div(B)=0, mais là au bord je touchait l'un des terme en le forçant à 0 (par exemple dBx/dx sur les bords x=cte) , l'autre (dBy/dy) n'avait aucune raison de l'etre et je n'y touchais pas.

Je me dis donc, qu'avec les conditions limites de mon post précédent, ou je ne touche a aucun terme de div(B), j'ai plus de chances de ne pas poser de probleme aux bords.

[FR119492]

Salut !

Il me semble bien évident que ce que je disais était que la fonction est périodique et de période la taille de la boite.

Ce n'est pas si évident que ça: j'ai rencontré des cas où la période était le double de la taille de la boîte.

Maintenant, pour clarifier les choses, est-ce que ta "boîte" est un objet matériel ( en cuivre, en fer, en verre ou en n'importe quoi d'autre) ou est-ce que c'est une frontière virtuelle, tes x allant en réalité de -infini à + infini ?

Jean-Marc Blanc

[Heimdall]

Salut,

Salut !
Ce n'est pas si évident que ça: j'ai rencontré des cas où la période était le double de la taille de la boîte.

Maintenant, pour clarifier les choses, est-ce que ta "boîte" est un objet matériel ( en cuivre, en fer, en verre ou en n'importe quoi d'autre) ou est-ce que c'est une frontière virtuelle, tes x allant en réalité de -infini à + infini ?

Hum, je suis pas sûr de savoir répondre à ça. Je ne simule pas une "boite" au sens physique du terme, mais mon domaine doit bien s'arreter quelque part. D'après les quelques articles que j'ai vu, on parle de "ideal conducting boundaries" qui impliquent en y=0 et y=Y, By=0, dBx/dy=dBz/dy=0

[Heimdall]

Re salut,

Bon j'ai essayé ce que je vous ai dit :

bord gauche et droit (x=0 ET x=X, y)
-----------------------------------------------
Bx(0,y) = Bx(X,y) périodique
By(0,y) = By(X,y) périodique
Bz(0,y) = Bz(X,y) périodique


Bord inférieur (x, y=0)
------------------------------
dBx/dy(x,0) = 0 (Neumann)
By(x,0) = 0 (Dirichlet)
dBz/dy(x,0) = 0 (Neumann)


bord supérieur (x, y=Y)
-------------------------------
dBx/dy(x,Y) = 0 (Neumann)
By(x,Y) = 0 (Dirichlet)
dBz/dy(x,Y) = 0 (Neumann)

Comment ? j'ai créé une grille buffer pour Sx, Sy, Sz qui :

- Suivant la direction x est entre x=0 et x=X
- suivant la direction y est entre y = 0 et y = 2Y

je fais :

Sy(x,y>Y) = -Sy(x,Y-y) (Sy(y) est impaire entre 0 et 2Y)
Sx(x,y>Y) = Sx(x,Y-y) (Sx(y) est paire entre 0 et 2Y)
Sz(x,y>Y) = Sz(x,Y-y) (Sz(y) est paire entre 0 et 2Y)

Puis j'effectue une transformée de fourier de ce nouveau S, et je résoud l'équation dans l'espace de fourier, c'est a dire en divisant la TF de S par (1+d^2+k^2).

Le prolongement impaire de Sy, permet d'assurer la condition de dirichlet (transformée en sinus) tandis que le prolongement de Sz et Sx permet d'assurer la condition de Neumann (transformée en cosinus)

La direction x restant inchangée entre x=0 et x=X, la transformée de fourier assure la condition périodique dans cette direction.

Et bien j'ai le même problème vous pouvez voir div(B ) ici http://nico.aunai.free.fr/divb.png

L'image est div(B ) représentée entre i = 0 et 64 (direction x), et j = 240 à 255 (direction y).

On voit clairement que div(B) augmente progressivement a mesure que l'on s'approche du bord y=Y. tandis qu'il n'y a aucun probleme dans la boite, et sur les cotés x=cte.

On pourrait croire que le probleme est encore assez faible, cependant l'erreur croit au court du temps et pose des problemes

[FR119492]

Salut !

Je crois que nous approchons du but, même si les informations que tu nous donnes sont parfaitement contradictoires: tu nous dis que les grandeurs sont périodiques en x mais tu donnes à X une valeur arbitraire. La période en x peut être n'importe quoi. Or les seules fonctions dont la période est indéterminée sont les constantes.

Alors, si tu veux qu'on puisse t'aider, énonce-nous ton problème d'une manière bien structurée:

1. Le problème tel qu'il t'a été posé initialement: domaine d'application, nature des phénomènes que tu étudies (ça semble être de l'électromagnétisme, donc équations de Maxwell); géométrie de l'objet étudié (forme et dimensions); matières constituant ce domaine (résistivité, permittivité, perméabilité); influences extérieures (tensions, courants, champs imposés) et leur évolution temporelles (régime constant, sinusoïdal, périodique, transitoire); et enfin ce que tu cherches.
2. Ce que tu as déjà essayé de faire, semble-t-il sans succès: hypothèses simplificatrices (3d ==> 2d; résistivité négligée; courant de déplacement négligé, etc.), tentatives infructueuses de résolution, etc.

Il est essentiel de bien séparer ce qui appartient au point 1, que nous devons accepter, et ce qui appartient au point 2, que nous pouvons rejeter.

Jean-Marc Blanc

[Heimdall]

bonjour,

Je crois que nous approchons du but, même si les informations que tu nous donnes sont parfaitement contradictoires: tu nous dis que les grandeurs sont périodiques en x mais tu donnes à X une valeur arbitraire. La période en x peut être n'importe quoi. Or les seules fonctions dont la période est indéterminée sont les constantes.

hum je comprends pas que tu aies tant de mal a accepter ce que je dis.
Quand on dit "conditions périodiques en x" cela signifie que *quelque soit* la taille de la boite dans la direction x, le bord x=X et x=0 seront fixés de telle maniere que les grandeurs sont périodiques de période X. C'est pas si dur à comprendre...

| i o xxxxxxxxxxxxxxxxx i o |
| i o xxxxxxxxxxxxxxxxx i o |
| i o xxxxxxxxxxxxxxxxx i o |
| i o xxxxxxxxxxxxxxxxx i o |
| i o xxxxxxxxxxxxxxxxx i o |


tu est bien d'accord ici que peut importe la valeur dans le domaine (symbolisé par les 'x') lorsque tu replie le domaine sur lui meme, les 'i' vont en face des 'i' de l'autre coté et les 'o' vont sur les 'o' de l'autre coté aussi. Le domaine est périodique, c'est tout !

Alors, si tu veux qu'on puisse t'aider, énonce-nous ton problème d'une manière bien structurée:

1. Le problème tel qu'il t'a été posé initialement: domaine d'application, nature des phénomènes que tu étudies (ça semble être de l'électromagnétisme, donc équations de Maxwell); géométrie de l'objet étudié (forme et dimensions); matières constituant ce domaine (résistivité, permittivité, perméabilité); influences extérieures (tensions, courants, champs imposés) et leur évolution temporelles (régime constant, sinusoïdal, périodique, transitoire); et enfin ce que tu cherches.
   
   

1. Bon alors je n'étais pas sûr de moi auparavant, mais maintenant si, il est complètement inutile de parler du code en général. Tout ce qu'il y a a savoir c'est comment inverser le systeme linéaire que j'ai donné ici, avec des conditions limites en y=0 et y=L qui respectent bien div(B)=0, et des conditions périodiques en x=0 et x=X.
   
   
   

2. Ce que tu as déjà essayé de faire, semble-t-il sans succès: hypothèses simplificatrices (3d ==> 2d; résistivité négligée; courant de déplacement négligé, etc.), tentatives infructueuses de résolution, etc.

[FR119492]

Alors, j'ai assez perdu de temps à essayer de t'aider. Maintenant, débrouille-toi sans moi.
Jean-Marc Blanc

[Heimdall]

Alors, j'ai assez perdu de temps à essayer de t'aider. Maintenant, débrouille-toi sans moi.
Jean-Marc Blanc

Dommage que tu le prenne comme ça, c'etait sympa d'essayer de m'aider tout de même. Peut-etre ai-je mal lu, mais tu n'arrete pas de dire que je dis n'importe quoi, même sur la définition de conditions périodiques. Je pense, meme si je n'ai apparement pas ton expérience dans le numérique, tout de meme savoir ce que je dis et ce que je veux dans le code. Je pense savoir ce qui est utile et pas utile.

et en particulier je pense de manière sûre qu'il nest pas utile de connaitre tout le code (je vais pas taper une doc de 500pages) pour savoir comment résoudre une équation bilinéaire sous condition d'avoir div(B).


Désolé de t'avoir dérangé, merci.

[Ehouarn]

Au vu du champs de div.B que tu obtiens, cela ressemble effectivement beaucoup à une incompatibilité entre conditions aux limites et div.B=0.

Pourtant à première vue, je ne vois pas d'incompatibilité flagrante au niveau de la méthode que tu emploies. Ceci dit, tu mentionnais le fait que tu résouts d'abord en S, puis que tu déduis ensuite B (via S=B-d^2laplacien(B)); est-ce que les conditions aux limites que tu prends (pour S) restent compatibles avec div B=0 ?

De même (c'est une question ouvert; je n'ai pas creusé), si div.B =0 est-ce-qu'on a alors div.laplacien(B)=0 ? Si c'est le cas, alors il pourrait être intéressant de regarder div.S pour cerner si ton problème est au niveau de la construction de S ou de la résolution de S=B-d^2laplacien(B).

[Heimdall]

Salut,

en fait, S est calculé par le rotationnel d'un champ, disons T. C'est fait par différences finies avec deux grilles. S est sur la grille "intérieure" (g1) et T sur une grille plus grande d'un dx/2 de chaque coté en x, et d'un dy/2 de chaque coté en y (appelons la g2), T est défini sur cette grande grille g2, de sorte que naturellement son rotationnel (et toute dérivée première) soit centrée sur la grille "intérieure" g1 ou se trouve S.

Il n'y a donc aucun besoin avec cette méthode, de donner des conditions "limites" sur S, étant donné que T, défini sur g2, permet de donner S partout sur la grille g1 sans restriction.

B est également défini sur cette grille intérieure. Auparavant il était lui aussi donné par le rotationnel d'un champ (E) défini sur g2, donc pas de probleme de conditions limites.

Maintenant, avec ma modification, ça n'est plus E mais T qui est sur g2, et du rotationnel de T je n'obtient pas B mais S, sur g1 donc. Il me faut passerpar la résolution du systeme linéaire pour déduire B, mais le prix a payer à cette intégration est de se donner des conditions limites pour B... Et bien sûr, celles ci doivent respecter div(B)=0.

J'ai déjà testé div(S), et ça donne 10^(-15) partout (i.e. sur g1), donc pas de probleme, je pense vraiment que le probleme vient de la résolution de S->B.

[Ehouarn]

D'accord pour l'utilisation judicieuse des grilles décalées. Mais pour le coup j'ai un petit doute sur l'utilisation d'une transformée en cosinus suivant y; les points de grilles associés ne doivent-ils pas normalement inclure la frontière (ce qui n'est pas le cas de la grille g1)?

Une autre interrogation en passant: tu utilises des différences finies pour calculer T et S, alors pourquoi ensuite passer par du Fourier pour ensuite calculer B; n'as-tu pas déjà tout l'arsenal pour la résolution du problème de Helmholtz via les différences finies sous la main?

[Heimdall]

D'accord pour l'utilisation judicieuse des grilles décalées. Mais pour le coup j'ai un petit doute sur l'utilisation d'une transformée en cosinus suivant y; les points de grilles associés ne doivent-ils pas normalement inclure la frontière (ce qui n'est pas le cas de la grille g1)?

je suis pas sûr de comprendre ta question. L'équation de helmoltz comporte un terme en B et un terme en d2B/dx2, le premier est naturellement sur g1 (B est défini sur g1), est le second l'est aussi. La dérivée premiere d'une quantité de g1 est sur g2, la dérivée seconde est sur g1 aussi mais un cran plus bas, il manque le tour.

La frontière "physique" du domaine est donnée par les points du tour de g1, les points du tour du g2 sont a l'extérieur est permettent de fixer les conditions limites des quantités définies sur g2. Par exemple :


0 1 2 3 4 5 6
0 + o + o + o +
1 + o + o
2 + o + o
3 + o +


Si + sont les points de g2, et o les points de g1. on peut dire que la croix en (1,0) (ligne,col) reçoit la meme composante tangentielle que celle en (1,2) et une composante normale oposée, assurant une composante normale nulle en (1,1) et la continuité de la composante tangentielle.

Cela n'est possible que pour une quantité définie sur g2.

Pour la transformée en sinus/cosinus, elles permettent d'assurer la condition de dirichlet et neumann respectivement (cf. Hockney et aussi NR) vu que l'une est la projection du signal sur une base de fonctions nulles aux bords et l'autre sur une base de fonctions a dérivée nulle au bord.

J'ai testé mon code sur un signal (2D mais invariant suivant y pour que je sache résoudre analytiquement l'équation), avec les conditions de Dirichlet et Neumann :

Neumann
Dirichlet

Ca a l'air de marcher, petit bémol pour la page Dirichlet, ou l'ont voit que numériquement ça colle bof pour 64 points, mais c'est parce que j'ai pris un cos() en signal d'entrée, la transition est abrupte pour tomber à 0, les signaux a fourieriser dans mon code ne sont pas si dramatiques (heureusement car 1024 points à 100 particules par cellule ça ferait mal)

je suis pas sûr d'avoir répondu a ta question ?

Une autre interrogation en passant: tu utilises des différences finies pour calculer T et S, alors pourquoi ensuite passer par du Fourier pour ensuite calculer B; n'as-tu pas déjà tout l'arsenal pour la résolution du problème de Helmholtz via les différences finies sous la main?

C'était ma première idée également. Mais pour intégrer ce genre d'équation en différences finies, il n'y a pas (a ma connaissance) 50 méthodes. Soit des méthodes de relaxations (bof rapides, et faut les coder), soit inverser la matrice de l'opérateur différentiel (1-d^2Laplacien) en incluant les conditions limites, mais pour une grille 1024^2 ça fait mal, soit les méthodes de fourier. Dans le hockney et dans NR ils parlent de méthode cyclic reduction et facr... je ne m'y suis pas attardé car les articles faisant ce que je cherche a faire parlent de fast fourier transform (et parlent des mêmes conditions périodique en x et non périodiques en y).

En plus par fourier, ça ne semble pas dramatiquement compliqué vu que l'opérateur est positif (1+d^2*k^2) aucune divergence n'est à craindre. Le seul truc qui ralentit c'est de symétriser le probleme de façon pair ou impaire suivant que l'on veut neumann ou dirichlet respectivement (cf mes liens web)

Je fais actuellement tourner un run homogene (toutes les quantités initiales sont constantes dans l'espace) et je vois toujours div(B) grimper sur les bords, pas de façon dramatique mais tout de même

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
________________ time step :    32820 ________________
 
divergence b     : 0.00251546691251212018
(z) impulsion    :  -0.1340
total energy     : 505.60311319274313746064
magnetic energy  : 206.2169
bulk p.energy    :   0.3719
p. energy        : 149.0143
e. energy        : 150.0000
electric fluct   : 0.07793286777597577264
cfl              :   0.0445
min # p. p.  cell : 71, (30,17)

alors que je suis parti de :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
 
divergence b     : 0.00001487900575514910
(z) impulsion    :  -0.0625
total energy     : 505.82254201658315651002
magnetic energy  : 206.2959
bulk p.energy    :   0.3823
p. energy        : 149.1443
e. energy        : 150.0000
electric fluct   : 0.06329316274032975587
cfl              :   0.0436
min # p. p. cell : 79, (58,37)

Les chiffres annoncés sont des maxima, j'ai regardé la carte 2D de div(B) et ce sont bien les bords non périodiques qui grimpent, les bords périodiques n'ont aucun probleme. On voit tout de meme que les énergies sont assez constantes, rien ne se passe, comme prévu. Mais au bout d'un moment lorsque div(B) va atteindre environ 0.1 des merdes vont arriver et le code va planter (en partie a cause du champ électrique qui va grimper en fleche et faire sortir ue particule à 12 années lumières de la boite)

Je me dis donc, (sous réserve d'une erreur de code mais j'y crois bof vu mes tests analytiques) que les conditions non périodiques que je prend ne respectent pas intrinsèquement div(B) = 0. Ma question était donc de savoir s'il existant de telles conditions limites, et la méthode qui permet de les imposer.