diff --git a/examples/check_deriv/deriv.input b/examples/check_deriv/deriv.input index 6bc4280..9ada05e 100644 --- a/examples/check_deriv/deriv.input +++ b/examples/check_deriv/deriv.input @@ -36,19 +36,29 @@ WEIGHT_CURV= 0.3000000000000000 , CASE_CURV=4 , PENFUN_CURV=2 , - CURV_ALPHA= 0.0000000000000000 , + CURV_ALPHA= 1.0000000000000000 , CURV_BETA= 2.0000000000000000 , CURV_GAMMA= 2.0000000000000000 , - CURV_SIGMA= 0.0000000000000000 , + CURV_SIGMA= 1.0000000000000000 , CURV_K0= 10.000000000000000 , curv_k1= 0.0000000000000000 , curv_k1len=0 , - WEIGHT_TORS= 0.4000000000000000 , + WEIGHT_STRAIGHT = 0.3, + CASE_STRAIGHT=4 , + PENFUN_STR=2 , + STR_ALPHA= 0.0000000000000000 , + STR_BETA= 2.0000000000000000 , + STR_GAMMA= 2.0000000000000000 , + STR_SIGMA= 1.0000000000000000 , + STR_K0= 3.000000000000000 , + str_k1= 0.0000000000000000 , + str_k1len=0 , + WEIGHT_TORS= 0.000004000000000000000 , CASE_TORS=1 , PENFUN_TORS=1 , TORS0= 0.0000000000000000 , - TORS_ALPHA= 1.0000000000000000 , - TORS_BETA= 2.0000000000000000 , + TORS_ALPHA= 1.00000000000000000 , + TORS_BETA= 2.000000000000000000 , TORS_GAMMA= 1.0000000000000000 , WEIGHT_NISSIN= 0.5000000000000000 , PENFUN_NISSIN=2 , diff --git a/examples/ellipse_spline/ellipse.focus b/examples/ellipse_spline/ellipse.focus new file mode 100644 index 0000000..c5a6acf --- /dev/null +++ b/examples/ellipse_spline/ellipse.focus @@ -0,0 +1,54 @@ + # Total number of coils + 4 + #----------------- 1 --------------------------- + #coil_type coil_symm coil_name + 5 2 Mod_001 + #Nseg current Ifree Length Lfree target_length + 128 1.000000000000000E+06 0 7.023857038530602E+00 1 7.000000000000000E+00 + #NCP + 35 +#vector T +-0.0941184909289301 -0.0627585813650836 -0.0313792041122836 0. 0.0313594161927663 0.0626824096405133 0.0939578445769685 0.125181792988915 0.1563581651782159 0.1874981472652705 0.2186185119074332 0.2497390511301777 0.2808795300648952 0.3120566449684125 0.3432814681347697 0.3745577763856298 0.40588150907107 0.4372414186349166 0.4686207958877165 0.5 0.5313594161927664 0.5626824096405133 0.5939578445769684 0.625181792988915 0.6563581651782158 0.6874981472652704 0.7186185119074332 0.7497390511301776 0.7808795300648952 0.8120566449684125 0.8432814681347697 0.8745577763856297 0.9058815090710699 0.9372414186349164 0.9686207958877164 1. 1.0313594161927664 1.0626824096405134 1.0939578445769684 +#Control points for coils ( x; y; z;) +4.040661646477116 4.063106738204567 4.045690028595451 3.9890717021253668 3.8954134991031255 3.7683002497575737 3.612606977203668 3.434315029250432 3.2402830008882355 3.037980845222405 2.835198151355947 2.6397394920357713 2.4591204893990377 2.300277672834113 2.16930362009899 2.0712168313579116 2.0097737630895907 1.9873286713621412 2.0047453809712583 2.0613637074413416 2.1550219104635837 2.2821351598091346 2.4378284323630375 2.616120380316276 2.8101524086784693 3.012454564344305 3.215237258210763 3.4106959175309366 3.591314920167672 3.7501577367325947 3.8811317894677146 3.9792185782088008 4.040661646477116 4.063106738204567 4.045690028595451 +0.5335499789954127 0.5190180212605333 0.4997391974872207 0.476452326727876 0.450050253821303 0.4215462239886883 0.3920354225094194 0.3626529827110528 0.3345300214716733 0.3087494959809587 0.2863038166063621 0.2680561274590347 0.2547069528161413 0.2467675450929382 0.2445408457521664 0.2481105696053506 0.2573385910020257 0.2718705487369057 0.2911493725102178 0.3144362432695632 0.3408383161761357 0.3693423460087507 0.3988531474880194 0.4282355872863858 0.456358548525765 0.4821390740164804 0.5045847533910772 0.5228324425384042 0.5361816171812976 0.5441210249045008 0.5463477242452721 0.542778000392089 0.5335499789954127 0.5190180212605333 0.4997391974872207 +-0.1761573604108771 0.0285906135067584 0.2333148297191023 0.4301593664131261 0.6115641241260152 0.7705537247670766 0.9010067802538252 0.9978949275335708 1.057481244120726 1.0774686040641839 1.0570905138140947 0.9971399721859984 0.8999354954052526 0.7692269472797766 0.6100466197424133 0.4285128859070993 0.2315948305160467 0.0268468565984102 -0.1778773596139325 -0.3747218963079568 -0.5561266540208465 -0.7151162546619069 -0.8455693101486563 -0.9424574574284008 -1.0020437740155572 -1.0220311339590142 -1.0016530437089255 -0.9417025020808282 -0.8444980253000844 -0.7137894771746077 -0.5546091496372443 -0.3730754158019304 -0.1761573604108771 0.0285906135067584 0.2333148297191023 +#----------------- 2 --------------------------- + #coil_type coil_symm coil_name + 5 2 Mod_002 + #Nseg current Ifree Length Lfree target_length + 128 1.000000000000000E+06 0 7.019510531933204E+00 1 7.000000000000000E+00 + #NCP + 35 +#vect +-0.0939656916711386 -0.0626151436995167 -0.0312894020840421 0. 0.0312470327064046 0.0624521056980779 0.0936220304388914 0.1247689998825705 0.1559087309834539 0.1870580534332631 0.2182323226218692 0.2494430653946117 0.2806962267426744 0.3119912833297375 0.3433213477863449 0.3746742334857137 0.4060343083288613 0.4373848563004831 0.4687105979159577 0.4999999999999998 0.5312470327064044 0.5624521056980777 0.5936220304388913 0.6247689998825703 0.6559087309834538 0.687058053433263 0.7182323226218691 0.7494430653946117 0.7806962267426744 0.8119912833297377 0.8433213477863449 0.8746742334857137 0.9060343083288614 0.9373848563004833 0.9687105979159579 1. 1.0312470327064045 1.062452105698078 1.0936220304388915 +#Control points for coils ( x; y; z;) +3.782220854888081 3.803902111047818 3.787906291555557 3.734836366888116 3.646722610174999 3.5269477692031552 3.380118320021826 3.2118862377043755 3.028728532527057 2.8376942070991573 2.6461298503298534 2.4613956215507726 2.290583038212049 2.140245099748302 2.0161481872345908 1.9230540688403128 1.8645392704023915 1.842858014242658 1.8588538337349174 1.91192375840236 2.000037515115477 2.1198123560873166 2.2666418052686494 2.434873887586099 2.6180315927634186 2.8090659181913145 3.000630274960623 3.1853645037397005 3.356177087078427 3.506515025542173 3.6306119380558837 3.7237060564501645 3.782220854888081 3.803902111047818 3.787906291555557 +1.5394681444391916 1.5417095204013247 1.529460906063889 1.503188533515847 1.463898749502908 1.4131004904212334 1.3527476508050122 1.285163406848953 1.212949602413353 1.138885164178981 1.065818024683164 0.9965551449841671 0.9337550133224252 0.8798265893655219 0.8368381916827832 0.8064393595830494 0.7897982756915847 0.7875568997294515 0.7998055140668866 0.8260778866149295 0.8653676706278678 0.9161659297095418 0.9765187693257656 1.0441030132818219 1.1163168177174236 1.190381255951794 1.2634483954476121 1.332711275146608 1.3955114068083525 1.449439830765254 1.4924282284479937 1.5228270605477254 1.5394681444391916 1.5417095204013247 1.529460906063889 +-0.1366558746147862 0.0684910320650332 0.2735474469580533 0.4706318039670198 0.6521620824122193 0.811149882909334 0.941473193328422 1.0381158900559453 1.0973632312691644 1.1169448443836805 1.0961196745638 1.0357005135520447 0.9380185750966424 0.806830869691396 0.6471748645095206 0.4651762496997887 0.2678167393699383 0.0626698326901195 -0.1423865822029007 -0.3394709392118679 -0.5210012176570655 -0.6799890181541819 -0.8103123285732695 -0.9069550253007927 -0.9662023665140107 -0.9857839796285285 -0.9649588098086462 -0.9045396487968929 -0.8068577103414897 -0.6756700049362434 -0.5160139997543686 -0.3340153849446363 -0.1366558746147862 0.0684910320650332 0.2735474469580533 +#----------------- 3 --------------------------- + #coil_type coil_symm coil_name + 5 2 Mod_003 + #Nseg current Ifree Length Lfree target_length + 128 1.000000000000000E+06 0 7.011869402917031E+00 1 7.000000000000000E+00 + #NCP + 35 +#vect +-0.0934578745603949 -0.062232800357122 -0.0310853532354074 0. 0.0310483611607238 0.0620905585077029 0.0931583763694472 0.1242796286400254 0.1554738784721072 0.1867494827042303 0.2181023967136856 0.2495168757530436 0.2809679162496244 0.3124250404499627 0.3438568614661525 0.3752357769380955 0.406542125439605 0.4377671996428781 0.4689146467645927 0.5000000000000001 0.531048361160724 0.562090558507703 0.5931583763694472 0.6242796286400255 0.6554738784721074 0.6867494827042304 0.7181023967136857 0.7495168757530437 0.7809679162496245 0.8124250404499629 0.8438568614661525 0.8752357769380954 0.9065421254396051 0.937767199642878 0.9689146467645926 1. 1.0310483611607237 1.062090558507703 1.0931583763694472 + #Fourier harmonics for coils ( xc; xs; yc; ys; zc; zs) +3.263747075912993 3.27558290033659 3.25473071595179 3.2019888561408636 3.1193928240497066 3.0101339512442533 2.8784305912046646 2.729359208827525 2.568654320580759 2.4024865131537574 2.2372273071193316 2.079209026682914 1.9344873761233643 1.8086141174017258 1.7064270112608901 1.6318640019892776 1.5878085211442825 1.5759726967206862 1.5968248811054837 1.6495667409164148 1.732162773007566 1.841421645813025 1.973125005852611 2.12219638822975 2.2829012764765175 2.4490690839035203 2.6143282899379474 2.772346570374359 2.917068220933912 3.042941479655549 3.145128585796386 3.219691595067999 3.263747075912993 3.27558290033659 3.25473071595179 +2.4677684212976225 2.487444582624472 2.4829369886027894 2.454415339474702 2.4029809041164665 2.3306223211867656 2.24013484957551 2.135007875921686 2.019286856526717 1.8974162530007503 1.7740698140786035 1.653974197740895 1.5417316511781314 1.4416472955936483 1.3575664463402797 1.2927273207949215 1.2496344727735849 1.2299583114467378 1.234465905468417 1.2629875545965068 1.3144219899547405 1.386780572884443 1.4772680444956978 1.5823950181495228 1.6981160375444901 1.8199866410704588 1.9433330799926056 2.0634286963303112 2.1756712428930767 2.2757555984775606 2.3598364477309284 2.4246755732762897 2.4677684212976225 2.487444582624472 2.4829369886027894 +-0.1382407046480632 0.0654694444512397 0.2691305874450401 0.4648918615456787 0.6452076876344496 0.8031357437168898 0.9326081124105098 1.0286637069815374 1.0876336391640753 1.1072745480346227 1.086847415688495 1.0271410664184415 0.9304407406577014 0.8004431780305289 0.6421206742136912 0.4615376331549941 0.2656243297066295 0.0619141806073269 -0.1417469623864733 -0.337508236487113 -0.5178240625758824 -0.6757521186583234 -0.8052244873519425 -0.9012800819229707 -0.9602500141055086 -0.9798909229760554 -0.959463790629929 -0.8997574413598752 -0.8030571155991347 -0.6730595529719617 -0.5147370491551255 -0.3341540080964281 -0.1382407046480632 0.0654694444512397 0.2691305874450401 +#----------------- 4 --------------------------- + #coil_type coil_symm coil_name + 5 2 Mod_004 + #Nseg current Ifree Length Lfree target_length + 128 1.000000000000000E+06 0 7.040968063122333E+00 1 7.000000000000000E+00 + #NCP + 35 +#vect +-0.0928170897138889 -0.061753175812694 -0.0308307874325665 0. 0.0308035383064877 0.061648494854175 0.0925969817739731 0.1236948354889956 0.1549646838402125 0.1864028214832725 0.2179801732118991 0.2496469664125492 0.2813402665424754 0.3129932585163177 0.344545031884861 0.3759495974708375 0.4071829102861112 0.4382468241873061 0.4691692125674335 0.5 0.5308035383064877 0.561648494854175 0.592596981773973 0.6236948354889956 0.6549646838402126 0.6864028214832725 0.7179801732118991 0.7496469664125494 0.7813402665424755 0.8129932585163178 0.8445450318848611 0.8759495974708376 0.9071829102861111 0.938246824187306 0.9691692125674335 1. 1.0308035383064877 1.061648494854175 1.0925969817739731 + #Fourier harmonics for coils ( xc; xs; yc; ys; zc; zs) +2.5090067416524717 2.5080339440151405 2.482086450239106 2.432169805903562 2.3602277205462907 2.2690575254008314 2.162190730619583 2.043748134013839 1.9182778634230009 1.7905827754048584 1.6655422990261513 1.5479334197140522 1.4422555622737843 1.3525641275841165 1.2823172406968175 1.2342402114304483 1.210212790095597 1.2111855877329292 1.237133081508965 1.287049725844506 1.3589918112017803 1.450162006347236 1.5570288011284894 1.6754713977342306 1.8009416683250687 1.9286367563432132 2.0536772327219186 2.171286112034018 2.276963969474286 2.366655404163953 2.4369022910512523 2.4849793203176223 2.5090067416524717 2.5080339440151405 2.482086450239106 +3.2404137081472077 3.2680212475574932 3.263205904187617 3.226155901587392 3.158322437525959 3.0623531449248085 2.941975180574969 2.801838809879038 2.647331842542329 2.4843731523835895 2.3191918560682767 2.1580981922842786 2.007252267873987 1.872436941644527 1.7588409581625675 1.6708583923066895 1.6119112204669344 1.5843036810566478 1.5891190244265254 1.626169027026749 1.694002491088184 1.7899717836893314 1.9103497480391762 2.0504861187351033 2.2049930860718128 2.367951776230554 2.533133072545866 2.694226736329864 2.8450726607401564 2.979887986969614 3.093483970451576 3.181466536307449 3.2404137081472077 3.2680212475574932 3.263205904187617 +-0.1748637053437298 0.0261337100004621 0.2271307854779068 0.4203495007153265 0.5983252951637416 0.7542079399642371 0.8820271580968926 0.9769132275354635 1.035268898666622 1.0548918953484883 1.0350477796820592 0.9764928072840477 0.8814466903056677 0.7535159524588462 0.5975692360275753 0.419566073956101 0.2263406410566901 0.0253432257124984 -0.1756538497649468 -0.3688725650023656 -0.5468483594507808 -0.7027310042512768 -0.8305502223839317 -0.9254362918225032 -0.9837919629536614 -1.0034149596355286 -0.9835708439690979 -0.9250158715710881 -0.8299697545927066 -0.7020390167458868 -0.5460923003146142 -0.3680891382431409 -0.1748637053437298 0.0261337100004621 0.2271307854779068 diff --git a/examples/ellipse_spline/ellipse.input b/examples/ellipse_spline/ellipse.input new file mode 100644 index 0000000..1d19791 --- /dev/null +++ b/examples/ellipse_spline/ellipse.input @@ -0,0 +1,84 @@ +&focusin + IsQuiet = -1 ! -2 verbose and including unconstrained cost functions; -1: verbose; 0: normal; 1: concise + IsSymmetric = 2 ! 0: no stellarator symmetry enforced; 1: plasma periodicity enforced; 2: coil periodicity enforced + + input_surf = '../rotating_ellipse/plasma.boundary' + case_surface = 0 ! 0: general VMEC-like format (Rbc, Rbs, Zbc, Zbs); 1: read axis for knots + knotsurf = 0.200D-00 ! minor plasma radius for knototrans, only valid for case surface = 1 + ellipticity = 0.000D+00 ! ellipticity of plasma for knototrans, only valid for case surface = 1 + Nteta = 64 ! poloidal number for discretizing the surface + Nzeta = 64 ! toroidal number for discretizing the surface + + case_init = 0 ! -1: read coils.ext file; 0: read ext.focus file; 1: initialize with circular coils; 2: initialize with dipoles + input_coils = 'ellipse.focus' + case_coils = 1 ! 0: using piecewise linear representation; (not ready); 1: using Fourier series representation + Ncoils = 4 ! number of coils; only valid when case_init = 1 + init_current = 1.000D+06 ! initial coil currents (Amper); only valid when case_init = 1 + init_radius = 0.500D+00 ! initial coil radius (meter); only valid when case_init = 1 + IsVaryCurrent = 0 ! 0: all the currents fixed; 1: currents can be changed; overwritten by ext.focus + IsVaryGeometry = 1 ! 0: all the geometries fixed; 1: geometries can be changed; overwritten by ext.focus + NFcoil = 4 ! number of Fourier harmonics representing the coils; overwritten by ext.focus + Nseg = 128 ! number of coil segments for discretizing; overwritten by ext.focus + + IsNormalize = 1 ! 0: do not normalize coil parameters; 1: normalize; I = I/I0, x = x/R0; I0 & R0 are quadrtic mean values. + IsNormWeight = 1 ! 0: do not normalize the weights; 1: normalize the weights + case_bnormal = 0 ! 0: keep raw Bn error; 1: Bn residue normalized to local |B| + case_length = 1 ! 1: quadratic format, converging the target length; 2: exponential format, as short as possible + weight_bnorm = 1.000D+02 ! weight for real space Bn errors + weight_bharm = 0.000D+00 ! weight for Bnm harmonic errors + weight_tflux = 0.000D+00 ! weight for toroidal flux error + target_tflux = 0.000D+00 ! target for the toroidal flux + weight_ttlen = 1.000D+01 ! weight for coil length error + target_length = 7.000D+00 ! target value (or for normalization) of the coils length, if zero, automatically set to initial actual length + weight_specw = 0.000D+00 ! weight for spectral condensation error + weight_cssep = 0.000D+00 ! weight for coil-surface separation constraint + weight_inorm = 1.000D+00 ! weight for normalization of current. Larger weight makes the derivatives more important. + weight_gnorm = 1.000D+00 ! weight for normalization of geometric coefficients. Larger weight makes the derivatives more important. + + case_optimize = 1 ! -2: check the 2nd derivatives (not ready); -1: check the 1st derivatives; 0: no optimizations performed; 1: optimizing using the gradient (DF and/or CG); + exit_tol = 1.000D-04 ! Exit the optimizer if the percent change in the cost function over the last 5 steps is below this threshold + + DF_maxiter = 0 ! maximum iterations allowed for using Differential Flow (DF) + DF_xtol = 1.000D-08 ! relative error for ODE solver + DF_tausta = 0.000D+00 ! starting value of τ. Usually 0.0 is a good idea + DF_tauend = 1.000D-00 ! ending value of τ. The larger value of τend − τsta, the more optimized + + CG_maxiter = 50 ! maximum iterations allowed for using Conjugate Gradient (CG) + CG_xtol = 1.000D-08 ! the stopping criteria of finding minimum; if |dχ2/dX| < CG xtol, exit the optimization + CG_wolfe_c1 = 0.1 ! c1 value in the strong wolfe condition for line search, (0.0, 0.5) + CG_wolfe_c2 = 0.9 ! c2 value in the strong wolfe condition for line search; 0 < c1 < c2 < 1 + + LM_maxiter = 0 ! maximum iterations allowed for using Levenberg-Marquard (LM) + LM_xtol = 1.000D-08 ! if the relative error between two consecutivec iterates is at most xtol, the optimization terminates + LM_ftol = 1.000D-08 ! if both the actual and predicted relative reductions in the sum of squares are at most ftol, the optimization terminates; + LM_factor = 100.0 ! the initial step bound, which is set to the product of factor and the euclidean norm of diag*x if nonzero + + case_postproc = 3 ! 0: no extra post-processing; 1: evaluate the current coils; 2: write SPEC file; 3: perform Poincare plots; 4: calculates |B| Fourier harmonics in Boozer coordinates + save_freq = 1 ! frequency for writing output files; should be positive + save_coils = 1 ! flag for indicating whether write example.focus and example.coils + save_harmonics = 0 ! flag for indicating whether write example.harmonics + save_filaments = 0 ! flag for indicating whether write .example.filaments.xxxxxx + + update_plasma = 0 ! if == 1, write new example.plasma file with updated Bn harmonics. + pp_phi = 0.000D+00 ! toroidal plane for poincare plots, cylindrical angle phi = pp_phi*Pi + pp_raxis = 3.000D+00 ! pp_raxis, pp_zaxis are initial guesses for magnetic axis at the specified toroidal angle + pp_zaxis = 0.000D+00 ! If both zero, FOCUS will take the geometric center as initial guess + pp_rmax = 0.000D+00 ! pp_rmax, pp_zmax are the upper bounds for performing fieldline tracing + pp_zmax = 0.000D+00 ! FOCUS will start follow fieldlines at interpolation between (pp_raxis, pp_zaxis) and (pp_rmax, pp_zmax) + pp_ns = 10 ! number of following fieldlines + pp_maxiter = 1000 ! number of periods for each fieldline following + pp_xtol = 1.000D-06 ! tolarence of ODE solver during fieldline fowllowing +/ +&mgrid +! mgrid file dimensions +Rmin = 0.0 +Rmax = 0.0 +Zmin = 0.0 +Zmax = 0.0 +Pmin = 0.0 +Pmax = 6.283 +! resolutions +NR = 101 +NZ = 101 +NP = 72 +/ diff --git a/sources/Makefile b/sources/Makefile index f7ceef9..3dbf4ce 100644 --- a/sources/Makefile +++ b/sources/Makefile @@ -2,10 +2,12 @@ ############################################################################################################ - ALLFILES= globals initial surface rdsurf rdknot rdcoils packdof bfield \ + + ALLFILES= globals initial surface rdsurf rdknot rdbooz rdcoils packdof bfield \ fdcheck datalloc solvers descent congrad lmalg \ bmnharm bnormal torflux length surfsep nissin torsion coilsep curvature \ - saving diagnos specinp poinplot boozer wtmgrid focus + splines straight_out saving diagnos specinp poinplot boozer wtmgrid focus + HFILES= $(ALLFILES:=.f90) # raw source files FFILES= $(ALLFILES:=_m.F90) # Fortran 90 files PFILES= $(ALLFILES:=.pdf) # documentations diff --git a/sources/bfield.f90 b/sources/bfield.f90 index a237d08..ea74e0d 100644 --- a/sources/bfield.f90 +++ b/sources/bfield.f90 @@ -28,7 +28,7 @@ subroutine bfield0(icoil, x, y, z, tBx, tBy, tBz) ! Be careful if coils have different resolutions. !------------------------------------------------------------------------------------------------------ use globals, only: dp, coil, surf, Ncoils, Nteta, Nzeta, cosnfp, sinnfp, & - zero, myid, ounit, Nfp, pi2, half, two, one, bsconstant, MPI_COMM_FOCUS + zero, myid, ounit, Nfp, pi2, half, two, one, bsconstant, MPI_COMM_FOCUS,coil_type_spline use mpi implicit none @@ -72,7 +72,7 @@ subroutine bfield0(icoil, x, y, z, tBx, tBy, tBz) Bx = zero; By = zero; Bz = zero select case (coil(icoil)%type) ! Fourier coils - case(1) + case(1,coil_type_spline) ! Biot-Savart law do kseg = 0, coil(icoil)%NS-1 dlx = xx - coil(icoil)%xx(kseg) @@ -141,7 +141,7 @@ subroutine bfield1(icoil, x, y, z, tBx, tBy, tBz, ND) ! Discretizing factor is includeed; coil(icoil)%dd(kseg) !------------------------------------------------------------------------------------------------------ use globals, only: dp, coil, DoF, surf, NFcoil, Ncoils, Nteta, Nzeta, & - zero, myid, ounit, Nfp, one, bsconstant, cosnfp, sinnfp, MPI_COMM_FOCUS + zero, myid, ounit, Nfp, one, bsconstant, cosnfp, sinnfp, MPI_COMM_FOCUS,coil_type_spline use mpi implicit none @@ -190,7 +190,7 @@ subroutine bfield1(icoil, x, y, z, tBx, tBy, tBz, ND) Bx = zero; By = zero; Bz = zero select case (coil(icoil)%type) - case(1) + case(1,coil_type_spline) ! Fourier coils NS = coil(icoil)%NS do kseg = 0, NS-1 diff --git a/sources/curvature.f90 b/sources/curvature.f90 index c7d14d0..7b52593 100644 --- a/sources/curvature.f90 +++ b/sources/curvature.f90 @@ -20,13 +20,13 @@ subroutine curvature(ideriv) use globals, only: dp, zero, half, pi2, machprec, ncpu, myid, ounit, MPI_COMM_FOCUS, & coil, DoF, Ncoils, Nfixgeo, Ndof, curv, t1K, t2K, weight_curv, FouCoil, & - mcurv, icurv, LM_fvec, LM_fjac + mcurv, icurv, LM_fvec, LM_fjac,coil_type_spline,Splines use mpi implicit none INTEGER, INTENT(in) :: ideriv - INTEGER :: astat, ierr, icoil, idof, ND, NF, ivec + INTEGER :: astat, ierr, icoil, idof, ND, NF,NCP, ivec REAL :: curvAdd !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! @@ -38,7 +38,7 @@ subroutine curvature(ideriv) if( ideriv >= 0 ) then do icoil = 1, Ncoils - if( coil(icoil)%type .ne. 1 ) exit ! only for Fourier + if( (coil(icoil)%type .ne. 1) .AND. (coil(icoil)%type .ne. coil_type_spline) ) exit if( coil(icoil)%Lc /= 0 ) then ! if geometry is free call CurvDeriv0(icoil,curvAdd) curv = curv + curvAdd @@ -62,17 +62,20 @@ subroutine curvature(ideriv) idof = 0 do icoil = 1, Ncoils - if(coil(icoil)%type .ne. 1) exit ! only for Fourier + if( (coil(icoil)%type .ne. 1) .AND. (coil(icoil)%type .ne. coil_type_spline) ) exit ND = DoF(icoil)%ND - NF = FouCoil(icoil)%NF + if (coil(icoil)%type == 1) NF = FouCoil(icoil)%NF + if (coil(icoil)%type == coil_type_spline) NCP = Splines(icoil)%NCP if ( coil(icoil)%Ic /= 0 ) then !if current is free; idof = idof + 1 endif if ( coil(icoil)%Lc /= 0 ) then !if geometry is free; - call CurvDeriv1( icoil, t1K(idof+1:idof+ND), ND, NF ) + + if (coil(icoil)%type == coil_type_spline) call CurvDeriv1( icoil, t1K(idof+1:idof+ND), ND, NCP ) + if (coil(icoil)%type == 1) call CurvDeriv1( icoil, t1K(idof+1:idof+ND), ND, NF ) if (mcurv > 0) then ! L-M format of targets LM_fjac(icurv+ivec, idof+1:idof+ND) = weight_curv * t1K(idof+1:idof+ND) ivec = ivec + 1 @@ -98,7 +101,7 @@ subroutine CurvDeriv0(icoil,curvRet) use globals, only: dp, zero, pi2, ncpu, astat, ierr, myid, ounit, coil, NFcoil, Nseg, Ncoils, & case_curv, curv_alpha, curv_k0, curv_k1, curv_beta, curv_gamma, curv_sigma, penfun_curv, & - curv_k1len, MPI_COMM_FOCUS + curv_k1len, MPI_COMM_FOCUS,coil_type_spline,Splines use mpi implicit none @@ -108,7 +111,8 @@ subroutine CurvDeriv0(icoil,curvRet) !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! INTEGER :: kseg, NS REAL :: hypc, curv_hold, k1_use - REAL :: curvv(0:coil(icoil)%NS), xt(0:coil(icoil)%NS), yt(0:coil(icoil)%NS), zt(0:coil(icoil)%NS), xa(0:coil(icoil)%NS), ya(0:coil(icoil)%NS), za(0:coil(icoil)%NS) + REAL :: curvv(0:coil(icoil)%NS), xt(0:coil(icoil)%NS), yt(0:coil(icoil)%NS), & + zt(0:coil(icoil)%NS), xa(0:coil(icoil)%NS), ya(0:coil(icoil)%NS), za(0:coil(icoil)%NS) NS = coil(icoil)%NS xt(0:coil(icoil)%NS) = coil(icoil)%xt(0:coil(icoil)%NS) @@ -150,6 +154,8 @@ subroutine CurvDeriv0(icoil,curvRet) curvv = sqrt( (za*yt-zt*ya)**2 + (xa*zt-xt*za)**2 + (ya*xt-yt*xa)**2 ) / (xt**2+yt**2+zt**2)**(1.5) coil(icoil)%maxcurv = maxval(curvv) + coil(icoil)%curvature = curvv + if ( curv_k1len .eq. 1 ) then k1_use = pi2/coil(icoil)%Lo else @@ -180,38 +186,43 @@ subroutine CurvDeriv0(icoil,curvRet) if ( curvv(kseg) > k1_use ) then curv_hold = curv_hold + curv_sigma*( ( curvv(kseg) - k1_use )**curv_gamma ) endif + endif curvRet = curvRet + curv_hold*sqrt(xt(kseg)**2+yt(kseg)**2+zt(kseg)**2) enddo call lenDeriv0( icoil, coil(icoil)%L ) - curvRet = pi2*curvRet/(NS*coil(icoil)%L) + if (coil(icoil)%type==1) curvRet = pi2*curvRet/(NS*coil(icoil)%L) + if (coil(icoil)%type==coil_type_spline) curvRet = 1.0*curvRet/(NS*coil(icoil)%L) + return end subroutine CurvDeriv0 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -subroutine CurvDeriv1(icoil, derivs, ND, NF) !Calculate all derivatives for a coil +subroutine CurvDeriv1(icoil, derivs, ND, NC) !Calculate all derivatives for a coil use globals, only: dp, zero, pi2, coil, DoF, myid, ounit, Ncoils, & case_curv, curv_alpha, curv_k0, curv_k1, curv_beta, curv_gamma, & - curv_sigma, penfun_curv, curv_k1len, FouCoil, MPI_COMM_FOCUS + curv_sigma, penfun_curv, curv_k1len, FouCoil, MPI_COMM_FOCUS,Splines,coil_type_spline use mpi implicit none - INTEGER, intent(in) :: icoil, ND , NF + INTEGER, intent(in) :: icoil, ND , NC !NC is actually NCP for splines and NF for FouCoil REAL , intent(out) :: derivs(1:1, 1:ND) !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! INTEGER :: kseg, astat, ierr, doff, nff, NS, n - REAL :: xt, yt, zt, xa, ya, za, f1, f2, curvHold, penCurv, leng, hypc, hyps, curv_deriv, & + REAL :: xt, yt, zt, xa, ya, za, f1, f2,ncc, nss, ncp, nsp, curvHold, penCurv,curvH,leng, hypc, hyps, curv_deriv, & k1_use, rtxrax, rtxray, rtxraz REAL :: d1L(1:1, 1:ND) REAL,allocatable :: dxtdDoF(:,:), dytdDoF(:,:), dztdDoF(:,:), dxadDoF(:,:), dyadDoF(:,:), dzadDoF(:,:) FATAL( CurvDeriv1, icoil .lt. 1 .or. icoil .gt. Ncoils, icoil not in right range ) + + derivs(1,1:ND) = 0.0 NS = coil(icoil)%NS SALLOCATE(dxtdDoF, (0:NS,1:ND), zero) @@ -220,28 +231,42 @@ subroutine CurvDeriv1(icoil, derivs, ND, NF) !Calculate all derivatives for a co SALLOCATE(dxadDoF, (0:NS,1:ND), zero) SALLOCATE(dyadDoF, (0:NS,1:ND), zero) SALLOCATE(dzadDoF, (0:NS,1:ND), zero) - do n = 1,NF - dxtdDof(0:NS,n+1) = -1*FouCoil(icoil)%smt(0:NS,n) * n - dxtdDof(0:NS,n+NF+1) = FouCoil(icoil)%cmt(0:NS,n) * n - dytdDof(0:NS,n+2*NF+2) = -1*FouCoil(icoil)%smt(0:NS,n) * n - dytdDof(0:NS,n+3*NF+2) = FouCoil(icoil)%cmt(0:NS,n) * n - dztdDof(0:NS,n+4*NF+3) = -1*FouCoil(icoil)%smt(0:NS,n) * n - dztdDof(0:NS,n+5*NF+3) = FouCoil(icoil)%cmt(0:NS,n) * n - - dxadDof(0:NS,n+1) = -1*FouCoil(icoil)%cmt(0:NS,n) * n*n - dxadDof(0:NS,n+NF+1) = -1*FouCoil(icoil)%smt(0:NS,n) * n*n - dyadDof(0:NS,n+2*NF+2) = -1*FouCoil(icoil)%cmt(0:NS,n) * n*n - dyadDof(0:NS,n+3*NF+2) = -1*FouCoil(icoil)%smt(0:NS,n) * n*n - dzadDof(0:NS,n+4*NF+3) = -1*FouCoil(icoil)%cmt(0:NS,n) * n*n - dzadDof(0:NS,n+5*NF+3) = -1*FouCoil(icoil)%smt(0:NS,n) * n*n - enddo - + + select case (coil(icoil)%type) + case(1) + do n = 1,NC + dxtdDof(0:NS,n+1) = -1*FouCoil(icoil)%smt(0:NS,n) * n + dxtdDof(0:NS,n+NC+1) = FouCoil(icoil)%cmt(0:NS,n) * n + dytdDof(0:NS,n+2*NC+2) = -1*FouCoil(icoil)%smt(0:NS,n) * n + dytdDof(0:NS,n+3*NC+2) = FouCoil(icoil)%cmt(0:NS,n) * n + dztdDof(0:NS,n+4*NC+3) = -1*FouCoil(icoil)%smt(0:NS,n) * n + dztdDof(0:NS,n+5*NC+3) = FouCoil(icoil)%cmt(0:NS,n) * n + + dxadDof(0:NS,n+1) = -1*FouCoil(icoil)%cmt(0:NS,n) * n*n + dxadDof(0:NS,n+NC+1) = -1*FouCoil(icoil)%smt(0:NS,n) * n*n + dyadDof(0:NS,n+2*NC+2) = -1*FouCoil(icoil)%cmt(0:NS,n) * n*n + dyadDof(0:NS,n+3*NC+2) = -1*FouCoil(icoil)%smt(0:NS,n) * n*n + dzadDof(0:NS,n+4*NC+3) = -1*FouCoil(icoil)%cmt(0:NS,n) * n*n + dzadDof(0:NS,n+5*NC+3) = -1*FouCoil(icoil)%smt(0:NS,n) * n*n + enddo + case(coil_type_spline) + dxtdDof(0:NS,0:ND) = Splines(icoil)%db_dt(0:NS,0:ND) + dytdDof(0:NS,0:ND) = Splines(icoil)%db_dt(0:NS,0:ND) + dztdDof(0:NS,0:ND) = Splines(icoil)%db_dt(0:NS,0:ND) + + dxadDof(0:NS,0:ND) = Splines(icoil)%db_dt_2(0:NS,0:ND) + dyadDof(0:NS,0:ND) = Splines(icoil)%db_dt_2(0:NS,0:ND) + dzadDof(0:NS,0:ND) = Splines(icoil)%db_dt_2(0:NS,0:ND) + case default + FATAL( CurvDeriv1, .true. , invalid coil_type option ) + end select + derivs = zero d1L = zero call lenDeriv0( icoil, coil(icoil)%L ) leng = coil(icoil)%L call lenDeriv1( icoil, d1L(1:1,1:ND), ND ) - + if ( case_curv .eq. 1 ) then curv_alpha = 0.0 curv_sigma = 1.0 @@ -264,11 +289,12 @@ subroutine CurvDeriv1(icoil, derivs, ND, NF) !Calculate all derivatives for a co rtxrax = yt*za - zt*ya rtxray = zt*xa - xt*za rtxraz = xt*ya - yt*xa - f1 = sqrt( (xt*ya-xa*yt)**2 + (xt*za-xa*zt)**2 + (yt*za-ya*zt)**2 ) - f2 = (xt**2+yt**2+zt**2)**1.5 + f1 = sqrt( (xt*ya-xa*yt)**2 + (xt*za-xa*zt)**2 + (yt*za-ya*zt)**2 ); + f2 = (xt**2+yt**2+zt**2)**(1.5); curvHold = f1/f2 penCurv = 0.0 curv_deriv = 0.0 + if ( curv_k1len .eq. 1 ) then k1_use = pi2/coil(icoil)%Lo else @@ -299,9 +325,12 @@ subroutine CurvDeriv1(icoil, derivs, ND, NF) !Calculate all derivatives for a co + ( rtxrax*(dytdDof(kseg,1:ND)*za - dztdDof(kseg,1:ND)*ya + yt*dzadDof(kseg,1:ND) - zt*dyadDof(kseg,1:ND)) & + rtxray*(dztdDof(kseg,1:ND)*xa - dxtdDof(kseg,1:ND)*za + zt*dxadDof(kseg,1:ND) - xt*dzadDof(kseg,1:ND)) & + rtxraz*(dxtdDof(kseg,1:ND)*ya - dytdDof(kseg,1:ND)*xa + xt*dyadDof(kseg,1:ND) - yt*dxadDof(kseg,1:ND)) )/(f1*f2) )*curv_deriv + enddo - derivs = pi2*derivs/(NS*leng) + if (coil(icoil)%type == 1) derivs = pi2*derivs/(NS*leng) + if (coil(icoil)%type == coil_type_spline )derivs = derivs/(NS*leng) + !if( case_curv == 4 ) derivs = derivs/leng return diff --git a/sources/datalloc.f90 b/sources/datalloc.f90 index 02b3261..d64fde8 100644 --- a/sources/datalloc.f90 +++ b/sources/datalloc.f90 @@ -13,7 +13,7 @@ subroutine AllocData(type) INTEGER, intent(in) :: type - INTEGER :: icoil, idof, ND, NF, icur, imag, isurf, NS, mm, iseg + INTEGER :: icoil, idof, ND, NF,NCP, icur, imag, isurf, NS, mm, iseg,i,iter_cp REAL :: xtmp, mtmp, tt isurf = plasma @@ -82,6 +82,8 @@ subroutine AllocData(type) SALLOCATE( coil(icoil)%zb, (0:coil(icoil)%NS), zero ) SALLOCATE( coil(icoil)%dl, (0:coil(icoil)%NS), zero ) SALLOCATE( coil(icoil)%dd, (0:coil(icoil)%NS), zero ) + SALLOCATE( coil(icoil)%curvature, (0:coil(icoil)%NS), zero ) + SALLOCATE( coil(icoil)%straight, (0:coil(icoil)%NS), zero ) coil(icoil)%dd = pi2 / NS ! discretizing factor; case(2) #ifdef dposition @@ -93,6 +95,58 @@ subroutine AllocData(type) case(3) DoF(icoil)%ND = coil(icoil)%Lc * 1 ! number of DoF for background Bt, Bz SALLOCATE(DoF(icoil)%xdof, (1:DoF(icoil)%ND), zero) + case(coil_type_spline) + ! get number of DoF for each coil and allocate arrays; + NS = coil(icoil)%NS + NCP = Splines(icoil)%NCP + ND = (3*NCP) ! total variables for geometry + DoF(icoil)%ND = coil(icoil)%Lc * ND !# of DoF for icoil; + SALLOCATE(DoF(icoil)%xdof, (1:DoF(icoil)%ND), zero) + SALLOCATE(DoF(icoil)%xof , (0:coil(icoil)%NS-1, 1:ND), zero) + SALLOCATE(DoF(icoil)%yof , (0:coil(icoil)%NS-1, 1:ND), zero) + SALLOCATE(DoF(icoil)%zof , (0:coil(icoil)%NS-1, 1:ND), zero) + ! allocate and calculate trignometric functions for re-use + SALLOCATE( Splines(icoil)%basis_0, (0:NS-1, 0:NCP+2), zero ) + SALLOCATE( Splines(icoil)%basis_1, (0:NS-1, 0:NCP+1), zero ) + SALLOCATE( Splines(icoil)%basis_2, (0:NS-1, 0:NCP), zero ) + SALLOCATE( Splines(icoil)%basis_3, (0:NS-1, 0:NCP-1), zero ) + SALLOCATE( Splines(icoil)%db_dt , (0:NS-1, 0:NCP-1), zero ) + SALLOCATE( Splines(icoil)%db_dt_2, (0:NS-1, 0:NCP-1), zero ) + + do i =0, coil(icoil)%NS-1 + Splines(icoil)%eval_points(i) = 1.0*(i)/(coil(icoil)%NS) + enddo + + ! the derivatives of dx/dv + + call eval_spline_basis(icoil) !Compute spline basis and derivatives using analytic derivation + call eval_spline_basis1(icoil) + call eval_spline_basis2(icoil) + !call check_eval_basis(icoil) !Check derivatives of basis functions numerically + call enforce_spline_periodicity(icoil) !Forces the derivates in respect to the first three control points to be equal to the last three + + DoF(icoil)%xof(0:coil(icoil)%NS-1, 1: NCP) = Splines(icoil)%basis_3(0:coil(icoil)%NS-1, 0: NCP-1) !x/xc + DoF(icoil)%yof(0:coil(icoil)%NS-1, NCP+1:2*NCP) = Splines(icoil)%basis_3(0:coil(icoil)%NS-1, 0: NCP-1) !y/yc + DoF(icoil)%zof(0:coil(icoil)%NS-1, 2*NCP+1:3*NCP) = Splines(icoil)%basis_3(0:coil(icoil)%NS-1, 0: NCP-1) !z/zc + + ! allocate xyz data + SALLOCATE( coil(icoil)%xx, (0:coil(icoil)%NS), zero ) + SALLOCATE( coil(icoil)%yy, (0:coil(icoil)%NS), zero ) + SALLOCATE( coil(icoil)%zz, (0:coil(icoil)%NS), zero ) + SALLOCATE( coil(icoil)%xt, (0:coil(icoil)%NS), zero ) + SALLOCATE( coil(icoil)%yt, (0:coil(icoil)%NS), zero ) + SALLOCATE( coil(icoil)%zt, (0:coil(icoil)%NS), zero ) + SALLOCATE( coil(icoil)%xa, (0:coil(icoil)%NS), zero ) + SALLOCATE( coil(icoil)%ya, (0:coil(icoil)%NS), zero ) + SALLOCATE( coil(icoil)%za, (0:coil(icoil)%NS), zero ) + SALLOCATE( coil(icoil)%dl, (0:coil(icoil)%NS), zero ) + SALLOCATE( coil(icoil)%dd, (0:coil(icoil)%NS), zero ) + SALLOCATE( coil(icoil)%curvature, (0:coil(icoil)%NS), zero ) + SALLOCATE( coil(icoil)%straight, (0:coil(icoil)%NS), zero ) + + coil(icoil)%dd = 1.0/(coil(icoil)%NS) + + case default FATAL(AllocData, .true., not supported coil types) end select @@ -102,8 +156,10 @@ subroutine AllocData(type) do icoil = 1, Ncoils Ndof = Ndof + coil(icoil)%Ic + DoF(icoil)%ND - if (allocated(FouCoil)) then + if (coil(icoil)%type==1) then Tdof = Tdof + 1 + 6*(FouCoil(icoil)%NF)+3 + else if (coil(icoil)%type==coil_type_spline) then + Tdof = Tdof + 1 + 3*(Splines(icoil)%NCP) else Tdof = Tdof + coil(icoil)%Ic + DoF(icoil)%ND end if @@ -118,7 +174,7 @@ subroutine AllocData(type) SALLOCATE( xdof, (1:Ndof), zero ) ! dof vector; SALLOCATE( dofnorm, (1:Ndof), one ) ! dof normalized value vector; - SALLOCATE( evolution, (1:Nouts+1, 0:11), zero ) !evolution array; + SALLOCATE( evolution, (1:Nouts+1, 0:12), zero ) !evolution array; SALLOCATE( coilspace, (1:Nouts+1, 1:Tdof), zero ) ! all the coil parameters; ! determine dofnorm @@ -127,7 +183,7 @@ subroutine AllocData(type) Inorm = zero ; Mnorm = zero icur = 0 ; imag = 0 ! icur for coil current count, imag for dipole count do icoil = 1, Ncoils - if(coil(icoil)%type == 1 .or. coil(icoil)%type == 3 ) then + if(coil(icoil)%type == 1 .or. coil(icoil)%type == 3 .or. coil(icoil)%type == coil_type_spline ) then ! Fourier representation or central currents Inorm = Inorm + coil(icoil)%I**2 icur = icur + 1 @@ -159,7 +215,7 @@ subroutine AllocData(type) idof = 0 do icoil = 1, Ncoils - if(coil(icoil)%type == 1) then ! Fourier representation + if(coil(icoil)%type == 1 .or. coil(icoil)%type == coil_type_spline ) then ! Fourier representation if(coil(icoil)%Ic /= 0) then dofnorm(idof+1) = Inorm idof = idof + 1 @@ -248,7 +304,7 @@ subroutine AllocData(type) FATAL( AllocData, Ndof < 1, INVALID Ndof value ) SALLOCATE( t1E, (1:Ndof), zero ) !SALLOCATE( deriv, (1:Ndof, 0:7), zero ) - SALLOCATE( deriv, (1:Ndof, 0:9), zero ) + SALLOCATE( deriv, (1:Ndof, 0:10), zero ) ! Bnorm related; if (weight_bnorm > sqrtmachprec .or. weight_bharm > sqrtmachprec) then @@ -284,6 +340,10 @@ subroutine AllocData(type) SALLOCATE( t1C, (1:Ndof), zero ) endif + if (weight_straight > sqrtmachprec) then + SALLOCATE( t1Str, (1:Ndof), zero ) + endif + ! cssep needed; if (weight_cssep > sqrtmachprec) then SALLOCATE( t1S, (1:Ndof), zero ) diff --git a/sources/diagnos.f90 b/sources/diagnos.f90 index 20f8958..d63bd01 100644 --- a/sources/diagnos.f90 +++ b/sources/diagnos.f90 @@ -5,11 +5,14 @@ SUBROUTINE diagnos ! DATE: 07/13/2017 ! diagonose the coil performance !------------------------------------------------------------------------------------------------------ + use globals use mpi implicit none - INTEGER :: icoil, itmp, NF, idof, i, j, isurf, cs, ip, is, Npc, coilInd0, coilInd1, j0, per0, l0, ss0 + + INTEGER :: icoil, itmp, NF,NCP, idof, i, j, isurf, cs, ip, is, Npc, coilInd0, coilInd1, j0, per0, l0, ss0 + LOGICAL :: lwbnorm, l_raw REAL :: MaxCurv, AvgLength, MinCCdist, MinCPdist, tmp_dist, ReDot, ImDot, dum, AvgCurv, AvgTors, & MinLambda, MaxS, torsRet @@ -37,6 +40,7 @@ SUBROUTINE diagnos write(ounit, 100) "Nissin complexity", nissin write(ounit, 100) "Coil-coil separation", ccsep write(ounit, 100) "Coil-surf separation", cssep + write(ounit, 100) "Straight-out", str endif 100 format (8X, ": ", A20, " = ", ES15.8) @@ -56,7 +60,10 @@ SUBROUTINE diagnos coilspace(iout, idof+1:idof+NF ) = FouCoil(icoil)%ys(1:NF) ; idof = idof + NF coilspace(iout, idof+1:idof+NF+1) = FouCoil(icoil)%zc(0:NF) ; idof = idof + NF +1 coilspace(iout, idof+1:idof+NF ) = FouCoil(icoil)%zs(1:NF) ; idof = idof + NF -!!$ case default + case (coil_type_spline) + NCP = Splines(icoil)%NCP + coilspace(iout, idof+1:idof+NCP*3) = Splines(icoil)%Cpoints(0:3*NCP-1) ; idof = idof + 3*NCP +!!$ case default !!$ FATAL(descent, .true., not supported coil types) end select enddo @@ -67,7 +74,7 @@ SUBROUTINE diagnos if ( (case_length == 1) .and. (sum(coil(1:Ncoils)%Lo) < sqrtmachprec) ) coil(1:Ncoils)%Lo = one call length(0) do icoil = 1, Ncoils - if(coil(icoil)%type .ne. 1) cycle ! only for Fourier + if(coil(icoil)%type .ne. 1 .AND. (coil(icoil)%type .ne. coil_type_spline)) cycle ! only for Fourier AvgLength = AvgLength + coil(icoil)%L enddo AvgLength = AvgLength / Ncoils @@ -76,7 +83,7 @@ SUBROUTINE diagnos !-------------------------------coil maximum curvature---------------------------------------------------- MaxCurv = zero do icoil = 1, Ncoils - if(coil(icoil)%type .ne. 1) cycle ! only for Fourier + if(coil(icoil)%type .ne. 1 .AND. (coil(icoil)%type .ne. coil_type_spline)) cycle ! only for Fourier call CurvDeriv0(icoil,dum) !dummy return if (coil(icoil)%maxcurv .ge. MaxCurv) then MaxCurv = coil(icoil)%maxcurv @@ -93,7 +100,7 @@ SUBROUTINE diagnos !-------------------------------average coil curvature------------------------------------------------------- AvgCurv = zero do icoil = 1, Ncoils - if(coil(icoil)%type .ne. 1) exit ! only for Fourier + if(coil(icoil)%type .ne. 1 .AND. (coil(icoil)%type .ne. coil_type_spline)) exit ! only for Fourier call avgcurvature(icoil) AvgCurv = AvgCurv + coil(icoil)%avgcurv #ifdef DEBUG @@ -143,7 +150,7 @@ SUBROUTINE diagnos coilInd0 = 0 coilInd1 = 1 do icoil = 1, Ncoils - if(coil(icoil)%type .ne. 1) cycle ! only for Fourier + if(coil(icoil)%type .ne. 1 .AND. (coil(icoil)%type .ne. coil_type_spline)) cycle if(Ncoils .eq. 1) exit ! if only one coil ! Data for the first coil SALLOCATE(Atmp, (1:3,0:coil(icoil)%NS-1), zero) @@ -184,7 +191,7 @@ SUBROUTINE diagnos endif do itmp = 1, Ncoils ! skip self and non-Fourier coils - if (itmp == icoil .or. coil(icoil)%type /= 1) cycle + if (itmp == icoil .or. ((coil(icoil)%type .ne. 1) .AND. (coil(icoil)%type .ne. coil_type_spline))) cycle SALLOCATE(Btmp, (1:3,0:coil(itmp)%NS-1), zero) ! check if the coil is stellarator symmetric !select case (coil(icoil)%symm) @@ -231,7 +238,7 @@ SUBROUTINE diagnos minCPdist = infmax do icoil = 1, Ncoils - if(coil(icoil)%type .ne. 1) cycle ! only for Fourier + if(coil(icoil)%type .ne. 1 .AND. (coil(icoil)%type .ne. coil_type_spline)) cycle SALLOCATE(Atmp, (1:3,0:coil(icoil)%NS-1), zero) SALLOCATE(Btmp, (1:3,1:(Nteta*Nzeta)), zero) @@ -335,14 +342,15 @@ END SUBROUTINE diagnos subroutine avgcurvature(icoil) - use globals, only: dp, zero, pi2, ncpu, astat, ierr, myid, ounit, coil, NFcoil, Nseg, Ncoils + + use globals, only: dp, zero, pi2, ncpu, astat, ierr, myid, ounit, coil, NFcoil, Nseg, Ncoils,Splines,coil_type_spline use mpi implicit none INTEGER, INTENT(in) :: icoil REAL,allocatable :: davgcurv(:) - + call length(0) SALLOCATE(davgcurv, (0:coil(icoil)%NS-1), zero) davgcurv = sqrt( (coil(icoil)%za*coil(icoil)%yt-coil(icoil)%zt*coil(icoil)%ya)**2 & @@ -350,7 +358,10 @@ subroutine avgcurvature(icoil) + (coil(icoil)%ya*coil(icoil)%xt-coil(icoil)%yt*coil(icoil)%xa)**2 )& / ((coil(icoil)%xt)**2+(coil(icoil)%yt)**2+(coil(icoil)%zt)**2)**(1.5) davgcurv = davgcurv*sqrt(coil(icoil)%xt**2+coil(icoil)%yt**2+coil(icoil)%zt**2) - coil(icoil)%avgcurv = pi2*sum(davgcurv)/size(davgcurv) + + if (coil(icoil)%type == coil_type_spline) coil(icoil)%avgcurv = 1.0*sum(davgcurv)/size(davgcurv) + if (coil(icoil)%type == 1)coil(icoil)%avgcurv = pi2*sum(davgcurv)/size(davgcurv) + coil(icoil)%avgcurv = coil(icoil)%avgcurv/coil(icoil)%L return diff --git a/sources/fdcheck.f90 b/sources/fdcheck.f90 index 17de6ba..5f269e7 100644 --- a/sources/fdcheck.f90 +++ b/sources/fdcheck.f90 @@ -60,6 +60,7 @@ SUBROUTINE fdcheck( ideriv ) do idof = 1, Ndof ! perturbation will be relative. small = xdof(idof) * psmall + if (small < machprec) small = psmall !backward pertubation; tmp_xdof = xdof tmp_xdof(idof) = tmp_xdof(idof) - half * small diff --git a/sources/focus.f90 b/sources/focus.f90 index 0b79fb4..7c757c7 100644 --- a/sources/focus.f90 +++ b/sources/focus.f90 @@ -37,7 +37,7 @@ PROGRAM focus use globals, only: dp, ncpu, myid, ounit, ierr, astat, eunit, case_surface, case_coils, case_optimize, & case_postproc, xdof, time_initialize, time_optimize, time_postproc, & - version, MPI_COMM_FOCUS + version, MPI_COMM_FOCUS,plasma_surf_boozer use mpi !to enable gfortran mpi_wtime bugs; 07/20/2017 implicit none @@ -57,7 +57,7 @@ PROGRAM focus select case( case_surface ) - case( 0 ) ; call surface ! general format (VMEC-like) plasma boundary; + case( 0,plasma_surf_boozer ) ; call surface ! general format (VMEC-like) plasma boundary; case( 1 ) ; call rdknot ! knototran-like plasma boundary; !case( 2 ) ; call readwout ! read vmec output for plasma boundary and Boozer coordinates; for future; diff --git a/sources/globals.f90 b/sources/globals.f90 index ce91a4d..08a8384 100644 --- a/sources/globals.f90 +++ b/sources/globals.f90 @@ -15,7 +15,9 @@ module globals !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! - CHARACTER(10), parameter :: version='v0.15.00' ! version number + + CHARACTER(10), parameter :: version='v0.16.00' ! version number + !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! @@ -81,6 +83,7 @@ module globals ! surface related INTEGER :: IsSymmetric = 0 INTEGER :: case_surface = 0 + INTEGER,parameter :: plasma_surf_boozer = 5 CHARACTER(100) :: input_surf = 'plasma.boundary' ! surface file REAL :: knotsurf = 0.200D-00 REAL :: ellipticity = 0.000D+00 @@ -100,12 +103,9 @@ module globals ! normalizations INTEGER :: IsNormalize = 1 INTEGER :: IsNormWeight = 1 - REAL :: weight_inorm = 1.000D+00 - REAL :: weight_gnorm = 1.000D+00 - REAL :: weight_mnorm = 1.000D+00 ! Normal field error - REAL :: weight_bnorm = 0.000D+00 INTEGER :: case_bnormal = 0 + REAL :: weight_bnorm = 1.000D+00 ! Bmn resonant harmonics REAL :: weight_bharm = 0.000D+00 INTEGER :: bharm_jsurf = 0 @@ -129,6 +129,18 @@ module globals REAL :: curv_k0 = 1.000D+01 REAL :: curv_k1 = 0.000D+00 INTEGER :: curv_k1len = 0 + ! straight function + REAL :: weight_straight = 0.000D+00 + REAL :: coeff_disp_straight = 0.50D+00 + INTEGER :: case_straight = 3 + INTEGER :: penfun_str = 1 + REAL :: str_alpha = 0.000D+00 + REAL :: str_beta = 2.000D+00 + REAL :: str_gamma = 2.000D+00 + REAL :: str_sigma = 0.000D+00 + REAL :: str_k0 = 1.000D+01 + REAL :: str_k1 = 0.000D+00 + INTEGER :: str_k1len = 0 ! coil torsion REAL :: weight_tors = 0.000D+00 INTEGER :: case_tors = 1 @@ -157,6 +169,12 @@ module globals REAL :: cssep_sigma = 0.000D+00 ! coil-coil separation REAL :: weight_ccsep = 0.000D+00 + REAL :: weight_inorm = 1.000D+00 + REAL :: weight_gnorm = 1.000D+00 + REAL :: weight_mnorm = 1.000D+00 + REAL :: origin_surface_x = 0.000D+00 + REAL :: origin_surface_y = 0.000D+00 + REAL :: origin_surface_z = 0.000D+00 INTEGER :: penfun_ccsep = 1 INTEGER :: ccsep_skip = 0 REAL :: r_delta = 1.000D-01 @@ -254,6 +272,20 @@ module globals curv_k0 ,& curv_k1 ,& curv_k1len ,& + weight_straight, & + coeff_disp_straight, & + case_straight , & + penfun_str ,& + str_alpha ,& + str_beta ,& + str_gamma ,& + str_sigma ,& + str_k0 ,& + str_k1 ,& + str_k1len ,& + origin_surface_x, & + origin_surface_y, & + origin_surface_z, & weight_tors ,& case_tors ,& penfun_tors ,& @@ -335,7 +367,7 @@ module globals !latex \subsection{surface and coils data} type toroidalsurface INTEGER :: Nteta, Nzeta, Nfou=0, Nfp=0, NBnf=0 - REAL , allocatable :: Rbc(:), Zbs(:), Rbs(:), Zbc(:), Bnc(:), Bns(:) + REAL , allocatable :: Rbc(:), Zbs(:), Rbs(:), Zbc(:), Bnc(:), Bns(:),Pmnc(:), Pmns(:) REAL , allocatable :: xx(:,:), yy(:,:), zz(:,:), nx(:,:), ny(:,:), nz(:,:), & xt(:,:), yt(:,:), zt(:,:), xp(:,:), yp(:,:), zp(:,:), & ds(:,:), bn(:,:), pb(:,:), & @@ -349,7 +381,7 @@ module globals REAL :: I=zero, L=zero, Lo, maxcurv, ox, oy, oz, mt, mp, Bt, Bz, avgcurv, & minlambda, maxs REAL , allocatable :: xx(:), yy(:), zz(:), xt(:), yt(:), zt(:), xa(:), ya(:), za(:), & - xb(:), yb(:), zb(:), dl(:), dd(:) + xb(:), yb(:), zb(:), dl(:), dd(:),curvature(:),straight(:) character(10) :: name end type arbitrarycoil @@ -358,6 +390,14 @@ module globals REAL , allocatable :: xc(:), xs(:), yc(:), ys(:), zc(:), zs(:), cmt(:,:), smt(:,:) end type FourierCoil + type SplineCoil + INTEGER :: NCP,NT,lwrk,ier + INTEGER,allocatable :: iwrk(:) + REAL :: FP + REAL , allocatable :: vect (:),basis_0(:,:),basis_1(:,:),basis_2(:,:),basis_3(:,:),Cpoints_fit(:),eval_points(:),wrk(:),weights(:), & + Data_points(:),Cpoints(:),db_dt(:,:),db_dt_2(:,:) + end type SplineCoil + type DegreeOfFreedom INTEGER :: ND REAL , allocatable :: xdof(:), xof(:,:), yof(:,:), zof(:,:) @@ -367,6 +407,7 @@ module globals type(toroidalsurface), target, allocatable :: surf(:) type(FourierCoil) , target, allocatable :: FouCoil(:) type(DegreeOfFreedom), target, allocatable :: DoF(:) + type(SplineCoil) , allocatable :: Splines(:) INTEGER :: Nfp = 1, symmetry = 0, surf_Nfp = 1 INTEGER :: plasma = 1, limiter = 1 @@ -384,8 +425,8 @@ module globals !latex \subsection{Optimization} ! General target functions; INTEGER :: iout, Nouts, LM_iter, LM_mfvec - INTEGER :: ibnorm = 0, ibharm = 0, itflux = 0, ittlen = 0, icssep = 0, icurv = 0, iccsep = 0, itors = 0, inissin = 0 ! starting number - INTEGER :: mbnorm = 0, mbharm = 0, mtflux = 0, mttlen = 0, mcssep = 0, mcurv = 0, mccsep = 0, mtors = 0, mnissin = 0 ! numbers of targets + INTEGER :: ibnorm = 0, ibharm = 0, itflux = 0, ittlen = 0, icssep = 0, icurv = 0, istr=0, iccsep = 0, itors = 0, inissin = 0 ! starting number + INTEGER :: mbnorm = 0, mbharm = 0, mtflux = 0, mttlen = 0, mcssep = 0, mcurv = 0, mstr=0, mccsep = 0, mtors = 0, mnissin = 0 ! numbers of targets REAL :: chi, discretefactor, sumDE REAL , allocatable :: t1E(:), t2E(:,:), evolution(:,:), coilspace(:,:), deriv(:,:) REAL , allocatable :: LM_fvec(:), LM_fjac(:,:) @@ -415,6 +456,8 @@ module globals ! nissin complexity REAL :: nissin REAL , allocatable :: t1N(:), t2N(:,:) + REAL :: str + REAL , allocatable :: t1Str(:), t2Str(:,:) ! Coil-surface spearation INTEGER :: psurf = 1 ! the prevent surface label; default 1 is the plasma boundary REAL :: cssep @@ -439,7 +482,7 @@ module globals !coilsI stores the true currents, coil%I stores scaled current;? INTEGER, allocatable :: coilseg(:) character(LEN=20), allocatable :: coilname(:) - + INTEGER,parameter :: coil_type_spline = 5 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! !latex \subsection{Time counting} @@ -449,7 +492,7 @@ module globals !latex \subsection{Miscellaneous} REAL :: tmpw_bnorm, tmpw_tflux ,tmpt_tflux, tmpw_ttlen, tmpw_specw, tmpw_ccsep, tmpw_bharm, & - tmpw_curv, tmpw_tors, tmpw_nissin + tmpw_curv, tmpw_tors, tmpw_nissin, tmpw_str REAL :: overlap = 0.0 !tmp weight for saving to restart file REAL, allocatable :: mincc(:,:), coil_importance(:) diff --git a/sources/initial.f90 b/sources/initial.f90 index 208ab5b..0202f7f 100644 --- a/sources/initial.f90 +++ b/sources/initial.f90 @@ -78,12 +78,15 @@ !latex with radius of \inputvar{init\_radius} and a central infinitely long current. !latex Dipole magnetizations anc the central current are all set to \inputvar{init\_current}. !latex \ei -!latex +!latex \item \inputvar{coil\_type\_spline = 5} \\ +!latex \textit{Specify integer assigned to cubic B-spline representation} \\ +!latex !latex \item \inputvar{case\_coils = 1} \\ !latex \textit{Specify representation used for the initial coils, seen in \link{rdcoils}} \\ !latex \bi \vspace{-5mm} !latex \item[0:] using piecewise linear representation; (not ready) !latex \item[1:] using Fourier series representation; +!latex \item[coil_type_spline:] using cubic B-spline representation; !latex \ei !latex !latex \item \inputvar{Ncoils = 0} \\ @@ -183,8 +186,14 @@ !latex !latex \item \inputvar{weight\_Gnorm = 1.0} \\ !latex \textit{additional factor for normalizing geometric variables; the larger, the more optimized for -!latex coil shapes; seen in \link{rdcoils}} -!latex +!latex coil shapes; seen in \link{rdcoils}} +!latex +!latex \item \inputvar{origin_surface_x = 0} \\ +!latex \textit{X coordinate of surface center used in computing straight cost function} +!latex \item \inputvar{origin_surface_y = 0} \\ +!latex \textit{Y coordinate of surface center used in computing straight cost function} +!latex \item \inputvar{origin_surface_z = 0} \\ +!latex \textit{Z coordinate of surface center used in computing straight cost function} !latex \par \begin{tikzpicture} \draw[dashed] (0,1) -- (10,1); \end{tikzpicture} !latex !latex \item \inputvar{case\_optimize = 1} \\ @@ -498,6 +507,10 @@ subroutine check_input write(ounit, 1000) 'case_surface', case_surface, 'Read axis information for expanding plasma boundary.' if (IsQuiet < 0) write(ounit, '(8X,": knotsurf = " ES12.5 & & " ; ellipticity = " ES12.5)') knotsurf, ellipticity + case (plasma_surf_boozer) + inquire( file=trim(input_surf), exist=exist ) + FATAL( initial, .not.exist, plasma boundary file not provided ) + write(ounit, 1000) 'case_surface', case_surface, 'Read Plasma boundary in Boozer coordinates.' case default FATAL( initial, .true., selected surface type is not supported ) end select @@ -641,6 +654,21 @@ subroutine check_input case default FATAL( initial, .true., selected penfun_nissin is not supported ) end select + + select case ( case_straight ) + case ( 1 ) + if (IsQuiet < 1) write(ounit, 1000) 'case_straight', case_straight, 'Linear format of straight-out coil penalty.' + case ( 2 ) + if (IsQuiet < 1) write(ounit, 1000) 'case_straight', case_straight, 'Quadratic format of straight-out coil penalty.' + case ( 3 ) + if (IsQuiet < 1) write(ounit, 1000) 'case_straight', case_straight, 'Penalty function of straight-out coil.' + case ( 4 ) + if (IsQuiet < 1) write(ounit, 1000) 'case_straight', case_straight, 'Linear and Penalty function.' + case default + FATAL( initial, .true., selected case_stright is not supported ) + end select + + FATAL( initial, weight_bnorm < zero, illegal ) FATAL( initial, weight_bharm < zero, illegal ) @@ -707,6 +735,7 @@ subroutine check_input !tmpw_specw = weight_specw tmpw_ccsep = weight_ccsep tmpw_curv = weight_curv + tmpw_str = weight_straight !tmpw_cssep = weight_cssep tmpw_tors = weight_tors tmpw_nissin = weight_nissin diff --git a/sources/length.f90 b/sources/length.f90 index a631a54..c5b8c59 100644 --- a/sources/length.f90 +++ b/sources/length.f90 @@ -62,7 +62,7 @@ subroutine length(ideriv) use globals, only: dp, zero, half, pi2, machprec, ncpu, myid, ounit, & coil, DoF, Ncoils, Nfixgeo, Ndof, ttlen, t1L, t2L, case_length, & ittlen, mttlen, LM_fvec, LM_fjac, weight_ttlen, length_delta, & - MPI_COMM_FOCUS + MPI_COMM_FOCUS,coil_type_spline use mpi implicit none @@ -78,7 +78,7 @@ subroutine length(ideriv) if( ideriv >= 0 ) then do icoil = 1, Ncoils !only care about unique coils; - if(coil(icoil)%type == 1) then ! only for Fourier + if((coil(icoil)%type == 1) .OR. (coil(icoil)%type == coil_type_spline)) then !if( myid.ne.modulo(icoil-1,ncpu) ) cycle ! parallelization loop; call LenDeriv0(icoil, coil(icoil)%L) !RlBCAST( coil(icoil)%L, 1, modulo(icoil-1,ncpu) ) !broadcast each coil's length @@ -125,7 +125,7 @@ subroutine length(ideriv) idof = idof +1 endif if ( coil(icoil)%Lc /= 0 ) then !if geometry is free; - if(coil(icoil)%type .eq. 1) then ! only for Fourier + if((coil(icoil)%type .eq. 1) .OR. (coil(icoil)%type == coil_type_spline)) then ! only for Fourier ! calculate normalization if (case_length == 1) then norm(icoil) = (coil(icoil)%L - coil(icoil)%Lo) / coil(icoil)%Lo**2 ! quadratic; @@ -174,7 +174,7 @@ end subroutine length subroutine LenDeriv0(icoil, length) - use globals, only: dp, zero, coil, myid, ounit, Ncoils, MPI_COMM_FOCUS + use globals, only: dp, zero, coil, myid, ounit, Ncoils, MPI_COMM_FOCUS,coil_type_spline use mpi implicit none diff --git a/sources/packdof.f90 b/sources/packdof.f90 index 7d0d053..5fb6576 100644 --- a/sources/packdof.f90 +++ b/sources/packdof.f90 @@ -25,7 +25,7 @@ SUBROUTINE packdof(lxdof) ! DATE: 2017/03/19 !--------------------------------------------------------------------------------------------- use globals, only : dp, zero, myid, ounit, MPI_COMM_FOCUS, & - & case_coils, Ncoils, coil, DoF, Ndof, DoFnorm + & case_coils, Ncoils, coil, DoF, Ndof, DoFnorm,coil_type_spline, Splines use mpi implicit none @@ -78,7 +78,20 @@ SUBROUTINE packdof(lxdof) idof = idof + 1 endif !--------------------------------------------------------------------------------------------- - case default + case(coil_type_spline) + + if(coil(icoil)%Ic /= 0) then + lxdof(idof+1) = coil(icoil)%I + idof = idof + 1 + endif + + ND = DoF(icoil)%ND + if(coil(icoil)%Lc /= 0) then + lxdof(idof+1:idof+ND) = DoF(icoil)%xdof(1:ND) + idof = idof + ND + endif +!------------------------------------------------------------- + case default FATAL(packdof01, .true., not supported coil types) end select @@ -102,7 +115,7 @@ SUBROUTINE unpacking(lxdof) ! DATE: 2017/04/03 !--------------------------------------------------------------------------------------------- use globals, only: dp, zero, myid, ounit, MPI_COMM_FOCUS, & - & case_coils, Ncoils, coil, DoF, Ndof, DoFnorm + & case_coils, Ncoils, coil, DoF, Ndof, DoFnorm,coil_type_spline, Splines use mpi implicit none @@ -151,6 +164,20 @@ SUBROUTINE unpacking(lxdof) DoF(icoil)%xdof(1) = lxdof(idof+1) * dofnorm(idof+1) idof = idof + 1 endif + !--------------------------------------------------------------------------------------------- + case(coil_type_spline) + + if(coil(icoil)%Ic /= 0) then + coil(icoil)%I = lxdof(idof+1) * dofnorm(idof+1) + idof = idof + 1 + endif + + ND = DoF(icoil)%ND + if(coil(icoil)%Lc /= 0) then + DoF(icoil)%xdof(1:ND) = lxdof(idof+1:idof+ND) * dofnorm(idof+1:idof+ND) + idof = idof + ND + endif + !--------------------------------------------------------------------------------------------- case default FATAL(unpacking01, .true., not supported coil types) @@ -176,11 +203,11 @@ SUBROUTINE packcoil ! pack coil representation variables into DoF (only geometries without currents); ! DATE: 2017/03/25 !--------------------------------------------------------------------------------------------- - use globals, only: dp, zero, myid, ounit, case_coils, Ncoils, coil, FouCoil, DoF, MPI_COMM_FOCUS + use globals, only: dp, zero, myid, ounit, case_coils, Ncoils, coil, FouCoil, DoF, MPI_COMM_FOCUS,coil_type_spline, Splines use mpi implicit none - INTEGER :: icoil, idof, NF, ierr, astat + INTEGER :: icoil, idof, NF, ierr, astat,NCP FATAL( packcoil01, .not. allocated(coil) , illegal ) ! FATAL( packcoil, .not. allocated(FouCoil), illegal ) @@ -229,6 +256,18 @@ SUBROUTINE packcoil DoF(icoil)%xdof(idof+1) = coil(icoil)%Bz; idof = idof + 1 endif FATAL( packcoil05 , idof .ne. DoF(icoil)%ND, counting error in packing ) +!--------------------------------------------------------------------------------------------- + case(coil_type_spline) + ! get number of DoF for each coil and allocate arrays; + NCP = Splines(icoil)%NCP + DoF(icoil)%xdof = zero + + !pack Spline coefficients; + idof = 0 + if(coil(icoil)%Lc /= 0) then + DoF(icoil)%xdof(idof+1 : idof+NCP*3) = Splines(icoil)%Cpoints(0:3*NCP-1); idof = idof + NCP*3 + endif + FATAL( packcoil03 , idof .ne. DoF(icoil)%ND, counting error in packing ) !--------------------------------------------------------------------------------------------- case default FATAL(packcoil06, .true., not supported coil types) @@ -245,11 +284,11 @@ SUBROUTINE unpackcoil ! pack coil representation variables into DoF (only geometries without currents); ! DATE: 2017/03/25 !--------------------------------------------------------------------------------------------- - use globals, only: dp, zero, myid, ounit, case_coils, Ncoils, coil, FouCoil, DoF, MPI_COMM_FOCUS + use globals, only: dp, zero, myid, ounit, case_coils, Ncoils, coil, FouCoil, DoF, MPI_COMM_FOCUS,coil_type_spline, Splines use mpi implicit none - INTEGER :: icoil, idof, NF, ierr, astat + INTEGER :: icoil, idof, NF, ierr, astat,NCP FATAL( unpackcoil01, .not. allocated(coil) , illegal ) ! FATAL( unpackcoil, .not. allocated(FouCoil), illegal ) @@ -297,6 +336,16 @@ SUBROUTINE unpackcoil coil(icoil)%Bz = DoF(icoil)%xdof(idof+1) ; idof = idof + 1 endif FATAL( unpackcoil05 , idof .ne. DoF(icoil)%ND, counting error in packing ) +!--------------------------------------------------------------------------------------------- + case(coil_type_spline) + ! get number of DoF for each coil and allocate arrays; + NCP = Splines(icoil)%NCP + idof = 0 + if (coil(icoil)%Lc /= 0) then + !unpack Spline coefficients; + Splines(icoil)%Cpoints(0:NCP*3-1) = DoF(icoil)%xdof(idof+1 : idof+3*NCP) ; idof = idof + 3*NCP + endif + FATAL( unpackcoil03 , idof .ne. DoF(icoil)%ND, counting error in packing ) !--------------------------------------------------------------------------------------------- case default diff --git a/sources/rdbooz.f90 b/sources/rdbooz.f90 new file mode 100644 index 0000000..7472250 --- /dev/null +++ b/sources/rdbooz.f90 @@ -0,0 +1,405 @@ +!title (boundary) ! Plasma boundary in Boozer coordinates + +!latex \briefly{A Fourier representation for the plasma boundary in Boozer coordinates.} + +!latex \calledby{\link{xfocus}} +!l tex \calls{\link{}} + +!latex \tableofcontents + +!latex \subsection{Overview} +!latex \bi +!latex \item[1.] The transformation between the Boozer and VMEC angles is calculated by the code +!latex \href{https://bitbucket.org/lazerson_princeton/stellopt/wiki/BOOZ_XFORM}{\blu{BOOZ\_XFORM}}. +!latex By reading the output of BOOZ\_XFORM, we can easily get the information we need: $Rmnc\_B$, +!latex $Zmns\_B$ and $Pmns\_B$. The three quantities are all decomposed in Boozer coordinates $(\t_B, \z_B)$. +!latex (Please note: here and the following we assume stellarator symmetry, $lasym = F$.) +!latex +!latex \item[2.] To inverse the transformation and map to cylindrical coordinates, we need to do +!latex \be +!latex R(\t_B, \z_B) & = & \sum_{m,n} Rmnc\_B_{mn} \cos(m\t_B - n\z_B) \nonumber \\ +!latex \phi & = & \z_B + p(\t_B, \z_B) \\ +!latex Z(\t_B, \z_B) & = & \sum_{m,n} Zmns\_B_{mn} \sin(m\t_B - n\z_B) \nonumber +!latex \ee +!latex Here, $\phi$ is the cylindrical angle and $p(\t_B, \z_B) = \sum_{mn} Pmns\_B \sin(m\t_B - n\z_B)$. +!latex +!latex \item[3.] Now, we can deriv the transformation to Cartesian coordinates, as +!latex \be +!latex x(\t_B, \z_B) & = & R(\t_B, \z_B) * \cos(\phi) \nonumber \\ +!latex y(\t_B, \z_B) & = & R(\t_B, \z_B) * \sin(\phi) \\ +!latex z(\t_B, \z_B) & = & Z(\t_B, \z_B) \nonumber +!latex \ee +!latex +!latex \ei +!latex +!latex \subsection{More details} +!latex For numerical applications in FOCUS, more quantities are computed. +!latex \bi +!latex \item[1.] Tangential derivatives (in the following, subscripts $B$ in $\t_B$, $\z_B$ and $Rmnc\_B$ etc. are neglected.) +!latex \be +!latex \pdv{x}{\t} & = & \pdv{R}{\t} \cos(\phi) - R \sin(\phi) \pdv{\phi}{\t} \nonumber \\ +!latex & = & \sum -m Rmnc \sin(m\t-n\z) \cos(\phi) - R \sin(\phi) [\sum \ m Pmns \cos(m\t-n\z)] \\ +!latex \pdv{x}{\z} & = & \pdv{R}{\z} \cos(\phi) - R \sin(\phi) \pdv{\phi}{\z} \nonumber \\ +!latex & = & \sum \, n Rmnc \sin(m\t-n\z) \cos(\phi) - R \sin(\phi) [1 + \sum -n Pmns \cos(m\t-n\z)] +!latex \ee +!latex +!latex \be +!latex \pdv{y}{\t} & = & \pdv{R}{\t} \sin(\phi) + R \cos(\phi) \pdv{\phi}{\t} \nonumber \\ +!latex & = & \sum -m Rmnc \sin(m\t-n\z) \sin(\phi) + R \cos(\phi) [\sum \ m Pmns \cos(m\t-n\z)] \\ +!latex \pdv{y}{\z} & = & \pdv{R}{\z} \sin(\phi) + R \cos(\phi) \pdv{\phi}{\z} \nonumber \\ +!latex & = & \sum \, n Rmnc \sin(m\t-n\z) \sin(\phi) + R \cos(\phi) [1 + \sum -n Pmns \cos(m\t-n\z)] +!latex \ee +!latex +!latex \be +!latex \pdv{z}{\t} & = & \sum \ m Zmns \cos(m\t-n\z) \\ +!latex \pdv{z}{\z} & = & \sum -n Zmns \cos(m\t-n\z) +!latex \ee +!latex +!latex \item[2.] The \emph{unit} normal vector of surface is defined as +!latex \be +!latex \vect{n} = \frac{\vect{x_{\z}} \times \vect{x_{\t}}}{|\vect{x_{\z}} \times \vect{x_{\t}}|} \ . +!latex \ee +!latex If poloidal and toroidal angle are both in counter-clockwise direction, the normal vector will point outwards. +!latex +!latex \item[3.] The Jacobian is $J = |\vect{x_{\z}} \times \vect{x_{\t}}|$. +!latex +!latex \ei + +!latex \subsection{Numerical implemention} +!latex When \inputvar{Itopoloy = 2}, \FOCUS will read \emph{plasma.boundary}. +!latex An example is showing as, +!latex { \begin{verbatim} +!latex #bmn bNfp nbf +!latex 4 2 0 +!latex #plasma boundary +!latex # n m Rbc Rbs Zbc Zbs Pmnc Pmns +!latex 0 0 3.00 0.0 0.0 0.00 0.0 0.1 +!latex 0 1 0.30 0.0 0.0 -0.30 0.0 0.1 +!latex 1 0 0.00 0.0 0.0 -0.06 0.0 0.1 +!latex 1 1 -0.06 0.0 0.0 -0.06 0.0 0.1 +!latex #Bn harmonics +!latex # n m bnc bns +!latex 0 0 0.0 0.0 +!latex \end{verbatim} +!latex } + +!latex $Rbc$, $Zbs$ and $Pmnc$ are all zero in stellarator symmetry. +!latex Number of $B_n$ coefficients, $nbf$, should also be zero, except for RMP coils. +!latex + +!latex \section{Boozer coordinates integration in FOCUS} \label{booz_xform} +!latex BOOZ\_XFORM transforms the VMEC coordinates into the straight field line coordinates introduced by Boozer. +!latex Although Boozer coordinates are normally used for representing magnetic fields, we are going to use Boozer angles ($\t_B, \z_B$) parameterizing the flux surface, on which magnetic field lines are lying. +!latex The surface in cylindrical coordinates is transformed into Boozer coordinates by +!latex \begin{equation} +!latex \left \{ +!latex \begin{array}{lll} +!latex \ds R(\t_B, \z_B) = \sum_{m,n} R^B_{mn} \cos(m\t_B - n\z_B) ; \\ +!latex \ds \phi(\t_B, \z_B) = \z_B +\sum_{m,n} P^B_{mn} \sin(m\t_B - n\z_B) ; \\ +!latex \ds Z(\t_B, \z_B) = \sum_{m,n} Z^B_{mn} \sin(m\t_B - n\z_B) \ . +!latex \end{array} +!latex \right . +!latex \end{equation} +!latex Here, stellarator symmetry is assumed. +!latex Fourier coefficients, $R^B_{mn}$, $Z^B_{mn}$ and $P^B_{mn}$, are provided in BOOZ\_XFORM outputs. +!latex +!latex FOCUS uses Cartesian coordinates $\vect{x}(\t_B, \z_B)$, therefore, the flux surface is parameterized as +!latex \begin{equation} +!latex \left \{ +!latex \begin{array}{lll} +!latex x(\t_B, \z_B) = R(\t_B, \z_B) \cos(\phi) ; \\ +!latex y(\t_B, \z_B) = R(\t_B, \z_B) \sin(\phi) ; \\ +!latex z(\t_B, \z_B) = Z(\t_B, \z_B) \ . +!latex \end{array} +!latex \right . +!latex \end{equation} +!latex Then the tangential derivatives are +!latex \begin{align} +!latex \pdv{x}{\t_B} = & \pdv{R}{\t_B} \cos(\phi) - R \sin(\phi) \pdv{\phi}{\t_B} ;\\ +!latex \pdv{x}{\z_B} = & \pdv{R}{\z_B} \cos(\phi) - R \sin(\phi) \pdv{\phi}{\z_B} . +!latex \end{align} +!latex Likewise, $\partial y / \partial {\t_B}$, $\partial y / \partial {\z_B}$, $\partial z / \partial {\t_B}$ and $\partial z / \partial {\z_B}$ are computed. +!latex The surface normal is $\vect{n} = \vect{x}_{\z} \times \vect{x}_{\t}$ and so is the Jacobian $J = \frac{\partial (x, y, z)}{\partial (\t_B, \z_B)} = |\vect{x}_{\t} \times \vect{x}_{\z}|$. + +!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! + +subroutine rdbooz(filename, index) + + use globals, only: dp, zero, half, one, pi2, myid, ncpu, ounit, runit, ext, IsSymmetric, & + Nfp, surf, cosnfp, sinnfp, MPI_COMM_FOCUS, surf_Nfp, symmetry, plasma, & + IsQuiet, Nteta, Nzeta, discretefactor, tflux_sign + + use mpi + + implicit none + + CHARACTER*100, INTENT(IN) :: filename + INTEGER, INTENT(IN) :: index + +!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! + + INTEGER :: iosta, astat, ierr, ii, jj, imn, Nfou, NBnf, ip + REAL :: RR(0:2), ZZ(0:2), Phi(0:2), szeta, czeta, & + xx(1:3), xt(1:3), xz(1:3), ds(1:3), teta, zeta, arg, carg, sarg, dd, tmp, dz + +!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! + + if (myid .eq. 0) then + open (runit, file=trim(filename), status='old', action='read') + read (runit, *) !empty line + read (runit, *) surf(index)%Nfou, surf(index)%Nfp, surf(index)%NBnf !read dimensions + endif + + !Broadcast the values + IlBCAST(surf(index)%Nfou, 1, 0) + IlBCAST(surf(index)%Nfp, 1, 0) + IlBCAST(surf(index)%NBnf, 1, 0) + + FATAL(rdbooz, surf(index)%Nfou .le. 0, invalid) + FATAL(rdbooz, surf(index)%Nfp .le. 0, invalid) + FATAL(rdbooz, surf(index)%NBnf < 0, invalid) + Nfou = surf(index)%Nfou + NBnf = surf(index)%NBnf + + ! Allocate arrays + SALLOCATE(surf(index)%bim, (1:Nfou), 0) + SALLOCATE(surf(index)%bin, (1:Nfou), 0) + SALLOCATE(surf(index)%Rbc, (1:Nfou), zero) + SALLOCATE(surf(index)%Rbs, (1:Nfou), zero) + SALLOCATE(surf(index)%Zbc, (1:Nfou), zero) + SALLOCATE(surf(index)%Zbs, (1:Nfou), zero) + SALLOCATE(surf(index)%Pmnc, (1:Nfou), zero) + SALLOCATE(surf(index)%Pmns, (1:Nfou), zero) + + if (myid .eq. 0) then + read (runit, *) !empty line + read (runit, *) !empty line + do imn = 1, Nfou + read (runit, *) surf(index)%bin(imn), surf(index)%bim(imn), surf(index)%Rbc(imn), surf(index)%Rbs(imn), & + surf(index)%Zbc(imn), surf(index)%Zbs(imn), surf(index)%Pmnc(imn), surf(index)%Pmns(imn) + enddo + endif + + IlBCAST(surf(index)%bim(1:Nfou), Nfou, 0) + IlBCAST(surf(index)%bin(1:Nfou), Nfou, 0) + + surf(index)%bin(1:Nfou) = surf(index)%bin(1:Nfou)*surf(index)%Nfp ! Disable periodicity + + RlBCAST(surf(index)%Rbc(1:Nfou), Nfou, 0) + RlBCAST(surf(index)%Rbs(1:Nfou), Nfou, 0) + RlBCAST(surf(index)%Zbc(1:Nfou), Nfou, 0) + RlBCAST(surf(index)%Zbs(1:Nfou), Nfou, 0) + RlBCAST(surf(index)%Pmnc(1:Nfou), Nfou, 0) + RlBCAST(surf(index)%Pmns(1:Nfou), Nfou, 0) + + if (surf(index)%NBnf .gt. 0) then !read Bn terms + SALLOCATE(surf(index)%Bnim, (1:NBnf), 0) + SALLOCATE(surf(index)%Bnin, (1:NBnf), 0) + SALLOCATE(surf(index)%Bnc, (1:NBnf), zero) + SALLOCATE(surf(index)%Bns, (1:NBnf), zero) + + if (myid .eq. 0) then + read (runit, *) !empty line + read (runit, *) !empty line + do imn = 1, NBnf + read (runit, *) surf(index)%Bnin(imn), surf(index)%Bnim(imn), surf(index)%Bnc(imn), surf(index)%Bns(imn) + enddo + endif + + IlBCAST(surf(index)%Bnim(1:NBnf), NBnf, 0) + IlBCAST(surf(index)%Bnin(1:NBnf), NBnf, 0) + + surf(index)%Bnin(1:NBnf) = surf(index)%Bnin(1:NBnf)*surf(index)%Nfp ! Disable periodicity + + RlBCAST(surf(index)%Bnc(1:NBnf), NBnf, 0) + RlBCAST(surf(index)%Bns(1:NBnf), NBnf, 0) + endif + + if (myid .eq. 0) close (runit, iostat=iosta) + + IlBCAST(iosta, 1, 0) + + FATAL(rdbooz, iosta .ne. 0, error closing surface) + +!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! + if (myid == 0 .and. IsQuiet <= 0) then + write (ounit, *) "-----------Reading surface in Boozer coordinates--------------" + write (ounit, '("surface : The surface ", A," will be discretized in "I6" X "I6" elements.")') trim(filename), Nteta, Nzeta + write (ounit, '(8X": Nfou = " I06 " ; Nfp = " I06 " ; NBnf = " I06 " ;" )') surf(index)%Nfou, surf(index)%Nfp, surf(index)%NBnf + endif + + if (myid == 0 .and. IsQuiet <= -2) then ! very detailed output; + write (ounit, '(" : " 10x " : bim ="10i13 )') surf(index)%bim(1:Nfou) + write (ounit, '(" : " 10x " : bin ="10i13 )') surf(index)%bin(1:Nfou) + write (ounit, '(" : " 10x " : Rbc ="10es13.5)') surf(index)%Rbc(1:Nfou) + write (ounit, '(" : " 10x " : Rbs ="10es13.5)') surf(index)%Rbs(1:Nfou) + write (ounit, '(" : " 10x " : Zbc ="10es13.5)') surf(index)%Zbc(1:Nfou) + write (ounit, '(" : " 10x " : Zbs ="10es13.5)') surf(index)%Zbs(1:Nfou) + write (ounit, '("surface : " 10x " : Pmnc ="10es13.5)') surf(index)%Pmnc(1:Nfou) + write (ounit, '("surface : " 10x " : Pmns ="10es13.5)') surf(index)%Pmns(1:Nfou) + + if (NBnf > 0) then + write (ounit, '(" : " 10x " : Bnim ="10i13 )') surf(index)%Bnim(1:NBnf) + write (ounit, '(" : " 10x " : Bnin ="10i13 )') surf(index)%Bnin(1:NBnf) + write (ounit, '(" : " 10x " : Bnc ="10es13.5)') surf(index)%Bnc(1:NBnf) + write (ounit, '(" : " 10x " : Bns ="10es13.5)') surf(index)%Bns(1:NBnf) + endif + endif + +!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! + + surf(index)%Nteta = Nteta + surf(index)%Nzeta = Nzeta + + if (index == plasma) then + Nfp = surf(plasma)%Nfp + surf_Nfp = Nfp ! local surface Nfp + select case (IsSymmetric) + case (0) + surf_Nfp = 1 ! reset Nfp to 1; + symmetry = 0 + case (1) ! plasma and coil periodicity enabled; + symmetry = 0 + case (2) ! stellarator symmetry enforced; + symmetry = 1 + end select + + SALLOCATE(cosnfp, (1:Nfp), zero) + SALLOCATE(sinnfp, (1:Nfp), zero) + do ip = 1, Nfp + cosnfp(ip) = cos((ip - 1)*pi2/Nfp) + sinnfp(ip) = sin((ip - 1)*pi2/Nfp) + enddo + ! discretefactor = discretefactor/Nfp + surf(index)%Nzeta = Nzeta*surf_Nfp*2**symmetry ! the total number from [0, 2pi] + discretefactor = (pi2/surf(plasma)%Nteta)*(pi2/surf(plasma)%Nzeta) + endif + + SALLOCATE(surf(index)%xx, (0:Nteta - 1, 0:Nzeta - 1), zero) !x coordinates; + SALLOCATE(surf(index)%yy, (0:Nteta - 1, 0:Nzeta - 1), zero) !y coordinates + SALLOCATE(surf(index)%zz, (0:Nteta - 1, 0:Nzeta - 1), zero) !z coordinates + SALLOCATE(surf(index)%nx, (0:Nteta - 1, 0:Nzeta - 1), zero) !unit nx; + SALLOCATE(surf(index)%ny, (0:Nteta - 1, 0:Nzeta - 1), zero) !unit ny; + SALLOCATE(surf(index)%nz, (0:Nteta - 1, 0:Nzeta - 1), zero) !unit nz; + SALLOCATE(surf(index)%ds, (0:Nteta - 1, 0:Nzeta - 1), zero) !jacobian; + SALLOCATE(surf(index)%xt, (0:Nteta - 1, 0:Nzeta - 1), zero) !dx/dtheta; + SALLOCATE(surf(index)%yt, (0:Nteta - 1, 0:Nzeta - 1), zero) !dy/dtheta; + SALLOCATE(surf(index)%zt, (0:Nteta - 1, 0:Nzeta - 1), zero) !dz/dtheta; + SALLOCATE(surf(index)%pb, (0:Nteta - 1, 0:Nzeta - 1), zero) !target Bn; + SALLOCATE(surf(index)%xp, (0:Nteta - 1, 0:Nzeta - 1), zero) !dx/dzeta; + SALLOCATE(surf(index)%yp, (0:Nteta - 1, 0:Nzeta - 1), zero) !dy/dzeta; + SALLOCATE(surf(index)%zp, (0:Nteta - 1, 0:Nzeta - 1), zero) !dz/dzeta; + + surf(index)%vol = zero ! volume enclosed by plasma boundary + surf(index)%area = zero ! surface area + +! The center point value was used to discretize grid + do ii = 0, Nteta - 1 + teta = (ii + half)*pi2/surf(index)%Nteta + do jj = 0, Nzeta - 1 + zeta = (jj + half)*pi2/surf(index)%Nzeta + RR(0:2) = zero; ZZ(0:2) = zero; Phi(0:2) = zero + do imn = 1, surf(index)%Nfou + arg = surf(index)%bim(imn)*teta - surf(index)%bin(imn)*zeta + carg = cos(arg); sarg = sin(arg) + + RR(0) = RR(0) + surf(index)%Rbc(imn)*carg + surf(index)%Rbs(imn)*sarg + ZZ(0) = ZZ(0) + surf(index)%Zbc(imn)*carg + surf(index)%Zbs(imn)*sarg + Phi(0) = Phi(0) + surf(index)%Pmnc(imn)*carg + surf(index)%Pmns(imn)*sarg + + RR(1) = RR(1) + (-surf(index)%Rbc(imn)*sarg + surf(index)%Rbs(imn)*carg)*surf(index)%bim(imn) + ZZ(1) = ZZ(1) + (-surf(index)%Zbc(imn)*sarg + surf(index)%Zbs(imn)*carg)*surf(index)%bim(imn) + Phi(1) = Phi(1) + (-surf(index)%Pmnc(imn)*sarg + surf(index)%Pmns(imn)*carg)*surf(index)%bim(imn) + + RR(2) = RR(2) - (-surf(index)%Rbc(imn)*sarg + surf(index)%Rbs(imn)*carg)*surf(index)%bin(imn) + ZZ(2) = ZZ(2) - (-surf(index)%Zbc(imn)*sarg + surf(index)%Zbs(imn)*carg)*surf(index)%bin(imn) + Phi(2) = Phi(2) - (-surf(index)%Pmnc(imn)*sarg + surf(index)%Pmns(imn)*carg)*surf(index)%bin(imn) + enddo ! end of do imn; 12/12/2018; + + Phi(0) = Phi(0) + zeta + !Phi(1) = Phi(1) + Phi(2) = Phi(2) + one + szeta = sin(Phi(0)) + czeta = cos(Phi(0)) + xx(1:3) = (/RR(0)*czeta, RR(0)*szeta, ZZ(0)/) + xt(1:3) = (/RR(1)*czeta - RR(0)*szeta*Phi(1), & + RR(1)*szeta + RR(0)*czeta*Phi(1), & + ZZ(1)/) + xz(1:3) = (/RR(2)*czeta - RR(0)*szeta*Phi(2), & + RR(2)*szeta + RR(0)*czeta*Phi(2), & + ZZ(2)/) + ds(1:3) = -(/xt(2)*xz(3) - xt(3)*xz(2), & + xt(3)*xz(1) - xt(1)*xz(3), & + xt(1)*xz(2) - xt(2)*xz(1)/) !careful with the negative sign; means counterclockwise; + dd = sqrt(sum(ds(1:3)*ds(1:3))) + ! x, y, z coordinates for the surface; + surf(index)%xx(ii, jj) = xx(1) + surf(index)%yy(ii, jj) = xx(2) + surf(index)%zz(ii, jj) = xx(3) + ! dx/dt, dy/dt, dz/dt (dt for d theta) + surf(index)%xt(ii, jj) = xt(1) + surf(index)%yt(ii, jj) = xt(2) + surf(index)%zt(ii, jj) = xt(3) + ! dx/dp, dy/dp, dz/dp (dp for d zeta(phi)) + surf(index)%xp(ii, jj) = xz(1) + surf(index)%yp(ii, jj) = xz(2) + surf(index)%zp(ii, jj) = xz(3) + ! surface normal vectors and ds for the jacobian; + surf(index)%nx(ii, jj) = ds(1)/dd + surf(index)%ny(ii, jj) = ds(2)/dd + surf(index)%nz(ii, jj) = ds(3)/dd + surf(index)%ds(ii, jj) = dd + ! using Gauss theorom; V = \int_S x \cdot n dt dz + surf(index)%vol = surf(index)%vol + surf(index)%xx(ii, jj)*ds(1) & + & + surf(index)%yy(ii, jj)*ds(2) + surf(index)%zz(ii, jj)*ds(3) + ! surface area + surf(index)%area = surf(index)%area + surf(index)%ds(ii, jj) + enddo ! end of do jj; 14 Apr 16; + enddo ! end of do ii; 14 Apr 16; + + ! print volume and area + surf(index)%vol = abs(surf(index)%vol)/3*(pi2/surf(index)%Nteta)*(pi2/surf(index)%Nzeta) + surf(index)%area = abs(surf(index)%area)*(pi2/surf(index)%Nteta)*(pi2/surf(index)%Nzeta) + if (index == plasma) then + surf(index)%vol = surf(index)%vol*surf_Nfp*2**symmetry + surf(index)%area = surf(index)%area*surf_Nfp*2**symmetry + endif + + if (myid == 0 .and. IsQuiet <= 0) then + write (ounit, '(8X": Enclosed total surface volume ="ES12.5" m^3 ; area ="ES12.5" m^2." )') & + surf(index)%vol, surf(index)%area + endif + + ! check theta direction for the plasma surface and determine the toroidal flux sign + if (index == plasma) then + dz = surf(plasma)%zz(1, 0) - surf(plasma)%zz(0, 0) + if (dz > 0) then + ! counter-clockwise + if (myid == 0) write (ounit, '(8X": The theta angle used is counter-clockwise.")') + tflux_sign = -1 + else + ! clockwise + if (myid == 0) write (ounit, '(8X": The theta angle used is clockwise.")') + tflux_sign = 1 + endif + endif + + !calculate target Bn with input harmonics; 05 Jan 17; + if (surf(index)%NBnf > 0) then + do jj = 0, Nzeta - 1 + zeta = (jj + half)*pi2/surf(index)%Nzeta + do ii = 0, Nteta - 1 + teta = (ii + half)*pi2/surf(index)%Nteta + do imn = 1, surf(index)%NBnf + arg = surf(index)%Bnim(imn)*teta - surf(index)%Bnin(imn)*zeta + surf(index)%pb(ii, jj) = surf(index)%pb(ii, jj) + surf(index)%Bnc(imn)*cos(arg) + surf(index)%Bns(imn)*sin(arg) + enddo + enddo + enddo + endif + +!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! + + return + +end subroutine rdbooz diff --git a/sources/rdcoils.f90 b/sources/rdcoils.f90 index 677ffe4..b189743 100644 --- a/sources/rdcoils.f90 +++ b/sources/rdcoils.f90 @@ -54,6 +54,19 @@ !latex . !latex . !latex . +!latex#------------coil_type_spline-------------------------------- +!latex #coil_type symm coil_name +!latex 5 0 Mod_001 +!latex #Nseg current Ifree Length Lfree target_length +!latex 128 9.844910899889484E+05 1 5.889288927667147E+00 1 1.000000000000000E+00 +!latex #NCP +!latex 13 +!latex #Vector of knots +!latex -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 +!latex #Control Points Coordinates for coils ( x; y; z) +!latex 0.58 0.4354 0.5464 0.242 0.456546 0.1414 0.465 0.23455 0.3636 0.75746 0.58 0.4354 0.5464 +!latex 0.4353 1.3453 2.356 0.3453 0.3453 0.23444 2.454 5.2453 2.5654 0.53453 0.4353 1.3453 2.356 +!latex -1.243 2.242 0.234 1.25 2.245 0.68876 1.45 4.4543 1.3445 6.3545 -1.243 2.242 0.234 !latex !latex \end{raw} !latex \ei @@ -198,13 +211,14 @@ subroutine rdcoils allocate( FouCoil(1:Ncoils) ) allocate( coil(1:Ncoils) ) allocate( DoF(1:Ncoils) ) + allocate( Splines(1:Ncoils) ) ! master CPU read the coils if( myid==0 ) then do icoil = 1, Ncoils read( runit,*) read( runit,*) read( runit,*) coil(icoil)%type, coil(icoil)%symm, coil(icoil)%name - FATAL( rdcoils04, coil(icoil)%type < 1 .or. coil(icoil)%type > 3, illegal ) + FATAL( rdcoils04, (coil(icoil)%type < 1 .or. coil(icoil)%type > 3) .and. (coil(icoil)%type .NE. coil_type_spline), illegal ) FATAL( rdcoils05, coil(icoil)%symm < 0 .or. coil(icoil)%symm > 2, illegal ) if(coil(icoil)%type == 1) then ! Fourier representation read( runit,*) @@ -239,6 +253,37 @@ subroutine rdcoils read( runit,*) read( runit,*) coil(icoil)%Ic, coil(icoil)%I, coil(icoil)%Lc, coil(icoil)%Bz coil(icoil)%symm = 0 ! automatic reset to 0; might not be necessary; 2020/01/17 + else if(coil(icoil)%type == coil_type_spline) then ! Spline representation + read( runit,*) + read( runit,*) coil(icoil)%NS, coil(icoil)%I, coil(icoil)%Ic, & + & coil(icoil)%L, coil(icoil)%Lc, coil(icoil)%Lo + FATAL( rdcoils06, coil(icoil)%NS < 0 , illegal ) + FATAL( rdcoils07, coil(icoil)%Ic < 0 .or. coil(icoil)%Ic > 1, illegal ) + FATAL( rdcoils08, coil(icoil)%Lc < 0 .or. coil(icoil)%Lc > 1, illegal ) + FATAL( rdcoils09, coil(icoil)%L < zero , illegal ) + FATAL( rdcoils10, coil(icoil)%Lo < zero , illegal ) + read( runit,*) + read( runit,*) Splines(icoil)%NCP + Splines(icoil)%NT = Splines(icoil)%NCP + 4 + FATAL( rdcoils12, Splines(icoil)%NCP < 0 , illegal ) + FATAL( rdcoils12_2, Splines(icoil)%NT < 0 , illegal ) + FATAL( rdcoils12_3, Splines(icoil)%NT .NE. Splines(icoil)%NCP +4, illegal ) + FATAL( rdcoils12_4, Splines(icoil)%NT < 10, illegal ) + SALLOCATE( Splines(icoil)%vect, (0:Splines(icoil)%NT-1), zero ) + SALLOCATE( Splines(icoil)%eval_points, (0:coil(icoil)%NS-1), zero ) + SALLOCATE( Splines(icoil)%Cpoints, (0:Splines(icoil)%NCP * 3 - 1 ), zero ) + read( runit,*) + read( runit,*) Splines(icoil)%vect(0:Splines(icoil)%NT-1) + Splines(icoil)%vect(Splines(icoil)%NT-3) = 1.0 + Splines(icoil)%vect(4) - Splines(icoil)%vect(3) + Splines(icoil)%vect(Splines(icoil)%NT-2) = 1.0 + Splines(icoil)%vect(5) - Splines(icoil)%vect(3) + Splines(icoil)%vect(Splines(icoil)%NT-1) = 1.0 + Splines(icoil)%vect(6) - Splines(icoil)%vect(3) + Splines(icoil)%vect(0) = Splines(icoil)%vect(Splines(icoil)%NT-7) - Splines(icoil)%vect(Splines(icoil)%NT-4) + Splines(icoil)%vect(1) = Splines(icoil)%vect(Splines(icoil)%NT-6) - Splines(icoil)%vect(Splines(icoil)%NT-4) + Splines(icoil)%vect(2) = Splines(icoil)%vect(Splines(icoil)%NT-5) - Splines(icoil)%vect(Splines(icoil)%NT-4) + read( runit,*) + read( runit,*) Splines(icoil)%Cpoints(0:Splines(icoil)%NCP-1) + read( runit,*) Splines(icoil)%Cpoints(Splines(icoil)%NCP:Splines(icoil)%NCP*2-1) + read( runit,*) Splines(icoil)%Cpoints(Splines(icoil)%NCP*2:Splines(icoil)%NCP*3-1) else STOP " wrong coil type in rdcoils" call MPI_ABORT(MPI_COMM_FOCUS, 1, ierr) @@ -295,6 +340,25 @@ subroutine rdcoils if(coil(icoil)%Ic == 0) Nfixcur = Nfixcur + 1 ! if(coil(icoil)%Lc == 0) Nfixgeo = Nfixgeo + 1 Nfixgeo = Nfixgeo + 1 ! always treat as a fixed geometry + else if (coil(icoil)%type == coil_type_spline) then + IlBCAST( coil(icoil)%NS , 1 , 0 ) + RlBCAST( coil(icoil)%I , 1 , 0 ) + IlBCAST( coil(icoil)%Ic , 1 , 0 ) + RlBCAST( coil(icoil)%L , 1 , 0 ) + IlBCAST( coil(icoil)%Lc , 1 , 0 ) + RlBCAST( coil(icoil)%Lo , 1 , 0 ) + IlBCAST( Splines(icoil)%NCP , 1 , 0 ) + IlBCAST( Splines(icoil)%NT , 1 , 0 ) + if (.not. allocated(Splines(icoil)%vect) ) then + SALLOCATE( Splines(icoil)%vect, (0:Splines(icoil)%NT-1), zero ) + SALLOCATE( Splines(icoil)%eval_points, (0:coil(icoil)%NS-1), zero ) + SALLOCATE( Splines(icoil)%Cpoints, (0:Splines(icoil)%NCP * 3 -1 ), zero ) + endif + RlBCAST( Splines(icoil)%vect(0:Splines(icoil)%NT -1) , Splines(icoil)%NT , 0 ) + RlBCAST( Splines(icoil)%eval_points(0:coil(icoil)%NS-1) , coil(icoil)%NS , 0 ) + RlBCAST( Splines(icoil)%Cpoints(0:Splines(icoil)%NCP*3 - 1) , Splines(icoil)%NCP*3 , 0 ) + if(coil(icoil)%Ic == 0) Nfixcur = Nfixcur + 1 + if(coil(icoil)%Lc == 0) Nfixgeo = Nfixgeo + 1 else STOP " wrong coil type in rdcoils" call MPI_ABORT(MPI_COMM_FOCUS, 1, ierr) @@ -474,13 +538,13 @@ subroutine discoil(ifirst) ! if ifirst = 1, it will update all the coils; otherwise, only update free coils; ! date: 20170314 !--------------------------------------------------------------------------------------------- - use globals, only: dp, zero, pi2, myid, ounit, coil, FouCoil, Ncoils, DoF, MPI_COMM_FOCUS + use globals, only: dp, zero, pi2, myid, ounit, coil, FouCoil, Ncoils, DoF, MPI_COMM_FOCUS,Splines,coil_type_spline use mpi implicit none INTEGER, intent(in) :: ifirst - INTEGER :: icoil, iseg, mm, NS, NF, ierr, astat, ip + INTEGER :: icoil, iseg, mm, NS, NCP, NF, ierr, astat, ip REAL :: tt !------------------------------------------------------------------------------------------- do icoil = 1, Ncoils @@ -540,6 +604,54 @@ subroutine discoil(ifirst) case(3) + case( coil_type_spline ) + ! reset to zero for all the coils; + coil(icoil)%xx = zero + coil(icoil)%yy = zero + coil(icoil)%zz = zero + coil(icoil)%xt = zero + coil(icoil)%yt = zero + coil(icoil)%zt = zero + coil(icoil)%xa = zero + coil(icoil)%ya = zero + coil(icoil)%za = zero + NS = coil(icoil)%NS + NCP = Splines(icoil)%NCP ! allias variable for simplicity; + + !-------------------------enforce periodicity---------------------------------------------- + Splines(icoil)%Cpoints(NCP-3) = Splines(icoil)%Cpoints(0) + Splines(icoil)%Cpoints(2*NCP-3) = Splines(icoil)%Cpoints(NCP) + Splines(icoil)%Cpoints(3*NCP-3) = Splines(icoil)%Cpoints(2*NCP) + + Splines(icoil)%Cpoints(NCP-2) = Splines(icoil)%Cpoints(1) + Splines(icoil)%Cpoints(2*NCP-2) = Splines(icoil)%Cpoints(NCP+1) + Splines(icoil)%Cpoints(3*NCP-2) = Splines(icoil)%Cpoints(2*NCP+1) + + Splines(icoil)%Cpoints(NCP-1) = Splines(icoil)%Cpoints(2) + Splines(icoil)%Cpoints(2*NCP-1) = Splines(icoil)%Cpoints(NCP+2) + Splines(icoil)%Cpoints(3*NCP-1) = Splines(icoil)%Cpoints(2*NCP+2) + !-------------------------calculate coil data------------------------------------------------- + do iseg=0,NS-1 + coil(icoil)%xx(iseg) = SUM (Splines(icoil)%Cpoints(0:NCP-1)*Splines(icoil)%basis_3(iseg,0:NCP-1)) + coil(icoil)%yy(iseg) = SUM (Splines(icoil)%Cpoints(NCP:2*NCP-1)*Splines(icoil)%basis_3(iseg,0:NCP-1)) + coil(icoil)%zz(iseg) = SUM (Splines(icoil)%Cpoints(2*NCP:3*NCP-1)*Splines(icoil)%basis_3(iseg,0:NCP-1)) + coil(icoil)%xt(iseg) = SUM (Splines(icoil)%Cpoints(0:NCP-1)*Splines(icoil)%db_dt(iseg,0:NCP-1)) + coil(icoil)%yt(iseg) = SUM (Splines(icoil)%Cpoints(NCP:2*NCP-1)*Splines(icoil)%db_dt(iseg,0:NCP-1)) + coil(icoil)%zt(iseg) = SUM (Splines(icoil)%Cpoints(2*NCP:3*NCP-1)*Splines(icoil)%db_dt(iseg,0:NCP-1)) + coil(icoil)%xa(iseg) = SUM (Splines(icoil)%Cpoints(0:NCP-1)*Splines(icoil)%db_dt_2(iseg,0:NCP-1)) + coil(icoil)%ya(iseg) = SUM (Splines(icoil)%Cpoints(NCP:2*NCP-1)*Splines(icoil)%db_dt_2(iseg,0:NCP-1)) + coil(icoil)%za(iseg) = SUM (Splines(icoil)%Cpoints(2*NCP:3*NCP-1)*Splines(icoil)%db_dt_2(iseg,0:NCP-1)) + enddo + + coil(icoil)%xx(NS) = coil(icoil)%xx(0) + coil(icoil)%yy(NS) = coil(icoil)%yy(0) + coil(icoil)%zz(NS) = coil(icoil)%zz(0) + coil(icoil)%xt(NS) = coil(icoil)%xt(0) + coil(icoil)%yt(NS) = coil(icoil)%yt(0) + coil(icoil)%zt(NS) = coil(icoil)%zt(0) + coil(icoil)%xa(NS) = coil(icoil)%xa(0) + coil(icoil)%ya(NS) = coil(icoil)%ya(0) + coil(icoil)%za(NS) = coil(icoil)%za(0) case default FATAL(discoil, .true., not supported coil types) end select diff --git a/sources/rdsurf.f90 b/sources/rdsurf.f90 index 0f2967c..0d113ca 100644 --- a/sources/rdsurf.f90 +++ b/sources/rdsurf.f90 @@ -1,5 +1,5 @@ -!title (boundary) ! The plasma boundary is read from file +!title (boundary) ! The plasma boundary is read from file !latex \briefly{A Fourier representation for the plasma boundary is read from file } @@ -15,12 +15,12 @@ !latex \ds Z &= \sum Z_{mn}^c \, \cos(m\t - n\z) + Z_{mn}^s \, \sin(m\t - n\z) \nonumber !latex \ee !latex Usually, if we imply stellarator symmetry, then $R_{mn}^s$ and $Z_{mn}^c$ would be zero. -!latex +!latex !latex The positive driection for poloidal angle $\t$ is \red{counterclockwise} and for toroidal angle is also !latex \red{counterclockwise} from the top view. The positive surface normal should be pointed outwards. !latex \subsection{Variables} -!latex The Fourier harmonics of the plasma boundary are reqired in \verb+plasma.boundary+, -!latex and the format of this file is as follows: +!latex The Fourier harmonics of the plasma boundary are reqired in \verb+plasma.boundary+, +!latex and the format of this file is as follows: !latex \begin{raw} !latex Nfou ! integer: number of Fourier harmonics for the plasma boundary; !latex Nfp ! integer: number of field periodicity; @@ -37,7 +37,7 @@ !latex Zbs(1:bmn) ! real : cylindrical Z sine harmonics; !latex Bns(1:nbf) ! real : B normal sin harmonics; !latex Bnc(1:nbf) ! real : B normal cos harmonics; \end{raw} -!latex Note that immediately after reading (and broadcasting) +!latex Note that immediately after reading (and broadcasting) !latex \verb+bin+, the field periodicity factor is included, i.e. \verb+bin = bin * Nfp+. !latex \subsection{Sample file} !latex Example of the plasma.boundary file: @@ -63,285 +63,285 @@ !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! subroutine fousurf(filename, index) - use globals, only : dp, zero, half, pi2, myid, ounit, runit, IsQuiet, IsSymmetric, & - Nteta, Nzeta, surf, discretefactor, Nfp, plasma, symmetry, & + use globals, only: dp, zero, half, pi2, myid, ounit, runit, IsQuiet, IsSymmetric, & + Nteta, Nzeta, surf, discretefactor, Nfp, plasma, symmetry, & tflux_sign, cosnfp, sinnfp, MPI_COMM_FOCUS, surf_Nfp - use mpi - implicit none + use mpi + implicit none + + CHARACTER*100, INTENT(IN) :: filename + INTEGER, INTENT(IN) :: index - CHARACTER*100, INTENT(IN) :: filename - INTEGER, INTENT(IN) :: index - !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! - - INTEGER :: iosta, astat, ierr, ii, jj, imn, Nfou, Nbnf, ip - REAL :: RR(0:2), ZZ(0:2), szeta, czeta, xx(1:3), xt(1:3), xz(1:3), ds(1:3), & - teta, zeta, arg, dd, dz - - ! read the header - if( myid == 0 ) then - open(runit, file=trim(filename), status='old', action='read') - read(runit,*) !empty line - read(runit,*) surf(index)%Nfou, surf(index)%Nfp, surf(index)%NBnf !read dimensions - endif - - !Broadcast the values - IlBCAST( surf(index)%Nfou , 1, 0 ) - IlBCAST( surf(index)%Nfp , 1, 0 ) - IlBCAST( surf(index)%NBnf , 1, 0 ) - FATAL( rdsurf, surf(index)%Nfou <= 0, invalid ) - FATAL( rdsurf, surf(index)%Nfp <= 0, invalid ) - FATAL( rdsurf, surf(index)%NBnf < 0, invalid ) - Nfou = surf(index)%Nfou - NBnf = surf(index)%NBnf - - !allocate arrays - SALLOCATE( surf(index)%bim, (1:Nfou), 0 ) - SALLOCATE( surf(index)%bin, (1:Nfou), 0 ) - SALLOCATE( surf(index)%Rbc, (1:Nfou), zero ) - SALLOCATE( surf(index)%Rbs, (1:Nfou), zero ) - SALLOCATE( surf(index)%Zbc, (1:Nfou), zero ) - SALLOCATE( surf(index)%Zbs, (1:Nfou), zero ) - - if( myid == 0 ) then - read(runit,*) !empty line - read(runit,*) !empty line - do imn = 1, surf(index)%Nfou - read(runit,*) surf(index)%bin(imn), surf(index)%bim(imn), surf(index)%Rbc(imn), & - & surf(index)%Rbs(imn), surf(index)%Zbc(imn), surf(index)%Zbs(imn) - enddo - endif - - IlBCAST( surf(index)%bim(1:Nfou), surf(index)%Nfou, 0 ) - IlBCAST( surf(index)%bin(1:Nfou), surf(index)%Nfou, 0 ) - - surf(index)%bin(1:Nfou) = surf(index)%bin(1:Nfou) * surf(index)%Nfp !The full plasma; - - RlBCAST( surf(index)%Rbc(1:Nfou), surf(index)%Nfou, 0 ) - RlBCAST( surf(index)%Rbs(1:Nfou), surf(index)%Nfou, 0 ) - RlBCAST( surf(index)%Zbc(1:Nfou), surf(index)%Nfou, 0 ) - RlBCAST( surf(index)%Zbs(1:Nfou), surf(index)%Nfou, 0 ) - - !read Bnormal ditributions - if( surf(index)%NBnf > 0) then - SALLOCATE( surf(index)%Bnim, (1:NBnf), 0 ) - SALLOCATE( surf(index)%Bnin, (1:NBnf), 0 ) - SALLOCATE( surf(index)%Bnc , (1:NBnf), zero ) - SALLOCATE( surf(index)%Bns , (1:NBnf), zero ) - - if( myid == 0 ) then - read(runit,*) !empty line - read(runit,*) !empty line - do imn = 1, surf(index)%NBnf - read(runit,*) surf(index)%Bnin(imn), surf(index)%Bnim(imn), surf(index)%Bnc(imn), surf(index)%Bns(imn) - enddo - endif - - IlBCAST( surf(index)%Bnim(1:NBnf), surf(index)%NBnf, 0 ) - IlBCAST( surf(index)%Bnin(1:NBnf), surf(index)%NBnf, 0 ) - - !if (IsSymmetric == 0) - surf(index)%Bnin(1:NBnf) = surf(index)%Bnin(1:NBnf) * surf(index)%Nfp ! periodicity; - ! This should be consistent with bnftran; Before fully constructed the stellarator symmetry, - ! it's turned off; - - RlBCAST( surf(index)%Bnc(1:NBnf) , surf(index)%NBnf, 0 ) - RlBCAST( surf(index)%Bns(1:NBnf) , surf(index)%NBnf, 0 ) - endif - - if( myid == 0 ) close(runit,iostat=iosta) - - IlBCAST( iosta, 1, 0 ) - - FATAL( surface, iosta.ne.0, error closing the surface ) - - !-------------output for check------------------------------------------------------------------------- - if( myid == 0 .and. IsQuiet <= 0) then - write(ounit, *) "-----------Reading surface-----------------------------------" - write(ounit, '("surface : The surface ", A," will be discretized in "I6" X "I6" elements.")') trim(filename), Nteta, Nzeta - write(ounit, '(8X": Nfou = " I06 " ; Nfp = " I06 " ; NBnf = " I06 " ;" )') surf(index)%Nfou, surf(index)%Nfp, surf(index)%NBnf - endif - - if( myid == 0 .and. IsQuiet <= -2) then ! very detailed output; - write(ounit,'(" : " 10x " : bim ="10i13 )') surf(index)%bim(1:Nfou) - write(ounit,'(" : " 10x " : bin ="10i13 )') surf(index)%bin(1:Nfou) - write(ounit,'(" : " 10x " : Rbc ="10es13.5)') surf(index)%Rbc(1:Nfou) - write(ounit,'(" : " 10x " : Rbs ="10es13.5)') surf(index)%Rbs(1:Nfou) - write(ounit,'(" : " 10x " : Zbc ="10es13.5)') surf(index)%Zbc(1:Nfou) - write(ounit,'(" : " 10x " : Zbs ="10es13.5)') surf(index)%Zbs(1:Nfou) - if(Nbnf > 0) then - write(ounit,'(" : " 10x " : Bnim ="10i13 )') surf(index)%Bnim(1:NBnf) - write(ounit,'(" : " 10x " : Bnin ="10i13 )') surf(index)%Bnin(1:NBnf) - write(ounit,'(" : " 10x " : Bnc ="10es13.5)') surf(index)%Bnc (1:NBnf) - write(ounit,'(" : " 10x " : Bns ="10es13.5)') surf(index)%Bns (1:NBnf) - endif - endif - - surf(index)%Nteta = Nteta - surf(index)%Nzeta = Nzeta - - if (index == plasma) then - Nfp = surf(plasma)%Nfp - surf_Nfp = Nfp ! local surface Nfp - select case (IsSymmetric) - case ( 0 ) - surf_Nfp = 1 ! reset Nfp to 1; - symmetry = 0 - case ( 1 ) ! plasma and coil periodicity enabled; - symmetry = 0 - case ( 2 ) ! stellarator symmetry enforced; - symmetry = 1 - end select - - SALLOCATE( cosnfp, (1:Nfp), zero ) - SALLOCATE( sinnfp, (1:Nfp), zero ) - do ip = 1, Nfp - cosnfp(ip) = cos((ip-1)*pi2/Nfp) - sinnfp(ip) = sin((ip-1)*pi2/Nfp) - enddo - ! discretefactor = discretefactor/Nfp - surf(index)%Nzeta = Nzeta * surf_Nfp * 2**symmetry ! the total number from [0, 2pi] - discretefactor = (pi2/surf(plasma)%Nteta) * (pi2/surf(plasma)%Nzeta) - endif - - SALLOCATE( surf(index)%xx, (0:Nteta-1,0:Nzeta-1), zero ) !x coordinates; - SALLOCATE( surf(index)%yy, (0:Nteta-1,0:Nzeta-1), zero ) !y coordinates - SALLOCATE( surf(index)%zz, (0:Nteta-1,0:Nzeta-1), zero ) !z coordinates - SALLOCATE( surf(index)%nx, (0:Nteta-1,0:Nzeta-1), zero ) !unit nx; - SALLOCATE( surf(index)%ny, (0:Nteta-1,0:Nzeta-1), zero ) !unit ny; - SALLOCATE( surf(index)%nz, (0:Nteta-1,0:Nzeta-1), zero ) !unit nz; - SALLOCATE( surf(index)%ds, (0:Nteta-1,0:Nzeta-1), zero ) !jacobian; - SALLOCATE( surf(index)%xt, (0:Nteta-1,0:Nzeta-1), zero ) !dx/dtheta; - SALLOCATE( surf(index)%yt, (0:Nteta-1,0:Nzeta-1), zero ) !dy/dtheta; - SALLOCATE( surf(index)%zt, (0:Nteta-1,0:Nzeta-1), zero ) !dz/dtheta; - SALLOCATE( surf(index)%pb, (0:Nteta-1,0:Nzeta-1), zero ) !target Bn; - SALLOCATE( surf(index)%xp, (0:Nteta-1,0:Nzeta-1), zero ) !dx/dzeta; - SALLOCATE( surf(index)%yp, (0:Nteta-1,0:Nzeta-1), zero ) !dy/dzeta; - SALLOCATE( surf(index)%zp, (0:Nteta-1,0:Nzeta-1), zero ) !dz/dzeta; - - surf(index)%vol = zero ! volume enclosed by plasma boundary - surf(index)%area = zero ! surface area - - ! The center point value was used to discretize grid; - do ii = 0, Nteta-1 - teta = ( ii + half ) * pi2 / surf(index)%Nteta - do jj = 0, Nzeta-1 - zeta = ( jj + half ) * pi2 / surf(index)%Nzeta - RR(0:2) = zero ; ZZ(0:2) = zero - do imn = 1, surf(index)%Nfou - arg = surf(index)%bim(imn) * teta - surf(index)%bin(imn) * zeta - RR(0) = RR(0) + surf(index)%Rbc(imn) * cos(arg) + surf(index)%Rbs(imn) * sin(arg) - ZZ(0) = ZZ(0) + surf(index)%Zbc(imn) * cos(arg) + surf(index)%Zbs(imn) * sin(arg) - RR(1) = RR(1) + ( - surf(index)%Rbc(imn) * sin(arg) + surf(index)%Rbs(imn) * cos(arg) ) * surf(index)%bim(imn) - ZZ(1) = ZZ(1) + ( - surf(index)%Zbc(imn) * sin(arg) + surf(index)%Zbs(imn) * cos(arg) ) * surf(index)%bim(imn) - RR(2) = RR(2) - ( - surf(index)%Rbc(imn) * sin(arg) + surf(index)%Rbs(imn) * cos(arg) ) * surf(index)%bin(imn) - ZZ(2) = ZZ(2) - ( - surf(index)%Zbc(imn) * sin(arg) + surf(index)%Zbs(imn) * cos(arg) ) * surf(index)%bin(imn) - enddo ! end of do imn; 30 Oct 15; - szeta = sin(zeta) - czeta = cos(zeta) - xx(1:3) = (/ RR(0) * czeta, RR(0) * szeta, ZZ(0) /) - xt(1:3) = (/ RR(1) * czeta, RR(1) * szeta, ZZ(1) /) - xz(1:3) = (/ RR(2) * czeta, RR(2) * szeta, ZZ(2) /) & - + (/ - RR(0) * szeta, RR(0) * czeta, zero /) - ! minus sign for theta counterclockwise direction; - ds(1:3) = -(/ xt(2) * xz(3) - xt(3) * xz(2), & - xt(3) * xz(1) - xt(1) * xz(3), & - xt(1) * xz(2) - xt(2) * xz(1) /) - dd = sqrt( sum( ds(1:3)*ds(1:3) ) ) - ! x, y, z coordinates for the surface; - surf(index)%xx(ii,jj) = xx(1) - surf(index)%yy(ii,jj) = xx(2) - surf(index)%zz(ii,jj) = xx(3) - ! dx/dt, dy/dt, dz/dt (dt for d theta) - surf(index)%xt(ii,jj) = xt(1) - surf(index)%yt(ii,jj) = xt(2) - surf(index)%zt(ii,jj) = xt(3) - ! dx/dp, dy/dp, dz/dp (dp for d zeta(phi)) - surf(index)%xp(ii,jj) = xz(1) - surf(index)%yp(ii,jj) = xz(2) - surf(index)%zp(ii,jj) = xz(3) - ! surface normal vectors and ds for the jacobian; - surf(index)%nx(ii,jj) = ds(1) / dd - surf(index)%ny(ii,jj) = ds(2) / dd - surf(index)%nz(ii,jj) = ds(3) / dd - surf(index)%ds(ii,jj) = dd - ! using Gauss theorom; V = \int_S x \cdot n dt dz - surf(index)%vol = surf(index)%vol + surf(index)%xx(ii,jj) * ds(1) & - & + surf(index)%yy(ii,jj) * ds(2) + surf(index)%zz(ii,jj) * ds(3) - ! surface area - surf(index)%area = surf(index)%area + surf(index)%ds(ii,jj) - enddo ! end of do jj; 14 Apr 16; - enddo ! end of do ii; 14 Apr 16; - - ! print volume and area - surf(index)%vol = abs(surf(index)%vol)/3 * (pi2/surf(index)%Nteta) * (pi2/surf(index)%Nzeta) - surf(index)%area = abs(surf(index)%area) * (pi2/surf(index)%Nteta) * (pi2/surf(index)%Nzeta) - if (index == plasma) then - surf(index)%vol = surf(index)%vol * surf_Nfp * 2**symmetry - surf(index)%area = surf(index)%area * surf_Nfp * 2**symmetry - endif - - if( myid == 0 .and. IsQuiet <= 0) then - write(ounit, '(8X": Enclosed total surface volume ="ES12.5" m^3 ; area ="ES12.5" m^2." )') & - surf(index)%vol, surf(index)%area - endif - - ! check theta direction for the plasma surface and determine the toroidal flux sign - if (index == plasma) then - dz = surf(plasma)%zz(1,0) - surf(plasma)%zz(0,0) - if (dz > 0) then - ! counter-clockwise - if( myid == 0) write(ounit, '(8X": The theta angle used is counter-clockwise.")') - tflux_sign = -1 - else - ! clockwise - if( myid == 0) write(ounit, '(8X": The theta angle used is clockwise.")') - tflux_sign = 1 - endif - endif - - !calculate target Bn with input harmonics; 05 Jan 17; - if(surf(index)%NBnf > 0) then - do jj = 0, Nzeta-1 - zeta = ( jj + half ) * pi2 / surf(index)%Nzeta - do ii = 0, Nteta-1 - teta = ( ii + half ) * pi2 / surf(index)%Nteta - do imn = 1, surf(index)%NBnf - arg = surf(index)%Bnim(imn) * teta - surf(index)%Bnin(imn) * zeta - surf(index)%pb(ii,jj) = surf(index)%pb(ii,jj) + surf(index)%Bnc(imn)*cos(arg) + surf(index)%Bns(imn)*sin(arg) - enddo - enddo - enddo - endif - - return - + + INTEGER :: iosta, astat, ierr, ii, jj, imn, Nfou, Nbnf, ip + REAL :: RR(0:2), ZZ(0:2), szeta, czeta, xx(1:3), xt(1:3), xz(1:3), ds(1:3), & + teta, zeta, arg, dd, dz + + ! read the header + if (myid == 0) then + open (runit, file=trim(filename), status='old', action='read') + read (runit, *) !empty line + read (runit, *) surf(index)%Nfou, surf(index)%Nfp, surf(index)%NBnf !read dimensions + endif + + !Broadcast the values + IlBCAST(surf(index)%Nfou, 1, 0) + IlBCAST(surf(index)%Nfp, 1, 0) + IlBCAST(surf(index)%NBnf, 1, 0) + FATAL(rdsurf, surf(index)%Nfou <= 0, invalid) + FATAL(rdsurf, surf(index)%Nfp <= 0, invalid) + FATAL(rdsurf, surf(index)%NBnf < 0, invalid) + Nfou = surf(index)%Nfou + NBnf = surf(index)%NBnf + + !allocate arrays + SALLOCATE(surf(index)%bim, (1:Nfou), 0) + SALLOCATE(surf(index)%bin, (1:Nfou), 0) + SALLOCATE(surf(index)%Rbc, (1:Nfou), zero) + SALLOCATE(surf(index)%Rbs, (1:Nfou), zero) + SALLOCATE(surf(index)%Zbc, (1:Nfou), zero) + SALLOCATE(surf(index)%Zbs, (1:Nfou), zero) + + if (myid == 0) then + read (runit, *) !empty line + read (runit, *) !empty line + do imn = 1, surf(index)%Nfou + read (runit, *) surf(index)%bin(imn), surf(index)%bim(imn), surf(index)%Rbc(imn), & + & surf(index)%Rbs(imn), surf(index)%Zbc(imn), surf(index)%Zbs(imn) + enddo + endif + + IlBCAST(surf(index)%bim(1:Nfou), surf(index)%Nfou, 0) + IlBCAST(surf(index)%bin(1:Nfou), surf(index)%Nfou, 0) + + surf(index)%bin(1:Nfou) = surf(index)%bin(1:Nfou)*surf(index)%Nfp !The full plasma; + + RlBCAST(surf(index)%Rbc(1:Nfou), surf(index)%Nfou, 0) + RlBCAST(surf(index)%Rbs(1:Nfou), surf(index)%Nfou, 0) + RlBCAST(surf(index)%Zbc(1:Nfou), surf(index)%Nfou, 0) + RlBCAST(surf(index)%Zbs(1:Nfou), surf(index)%Nfou, 0) + + !read Bnormal ditributions + if (surf(index)%NBnf > 0) then + SALLOCATE(surf(index)%Bnim, (1:NBnf), 0) + SALLOCATE(surf(index)%Bnin, (1:NBnf), 0) + SALLOCATE(surf(index)%Bnc, (1:NBnf), zero) + SALLOCATE(surf(index)%Bns, (1:NBnf), zero) + + if (myid == 0) then + read (runit, *) !empty line + read (runit, *) !empty line + do imn = 1, surf(index)%NBnf + read (runit, *) surf(index)%Bnin(imn), surf(index)%Bnim(imn), surf(index)%Bnc(imn), surf(index)%Bns(imn) + enddo + endif + + IlBCAST(surf(index)%Bnim(1:NBnf), surf(index)%NBnf, 0) + IlBCAST(surf(index)%Bnin(1:NBnf), surf(index)%NBnf, 0) + + !if (IsSymmetric == 0) + surf(index)%Bnin(1:NBnf) = surf(index)%Bnin(1:NBnf)*surf(index)%Nfp ! periodicity; + ! This should be consistent with bnftran; Before fully constructed the stellarator symmetry, + ! it's turned off; + + RlBCAST(surf(index)%Bnc(1:NBnf), surf(index)%NBnf, 0) + RlBCAST(surf(index)%Bns(1:NBnf), surf(index)%NBnf, 0) + endif + + if (myid == 0) close (runit, iostat=iosta) + + IlBCAST(iosta, 1, 0) + + FATAL(surface, iosta .ne. 0, error closing the surface) + + !-------------output for check------------------------------------------------------------------------- + if (myid == 0 .and. IsQuiet <= 0) then + write (ounit, *) "-----------Reading surface-----------------------------------" + write (ounit, '("surface : The surface ", A," will be discretized in "I6" X "I6" elements.")') trim(filename), Nteta, Nzeta + write (ounit, '(8X": Nfou = " I06 " ; Nfp = " I06 " ; NBnf = " I06 " ;" )') surf(index)%Nfou, surf(index)%Nfp, surf(index)%NBnf + endif + + if (myid == 0 .and. IsQuiet <= -2) then ! very detailed output; + write (ounit, '(" : " 10x " : bim ="10i13 )') surf(index)%bim(1:Nfou) + write (ounit, '(" : " 10x " : bin ="10i13 )') surf(index)%bin(1:Nfou) + write (ounit, '(" : " 10x " : Rbc ="10es13.5)') surf(index)%Rbc(1:Nfou) + write (ounit, '(" : " 10x " : Rbs ="10es13.5)') surf(index)%Rbs(1:Nfou) + write (ounit, '(" : " 10x " : Zbc ="10es13.5)') surf(index)%Zbc(1:Nfou) + write (ounit, '(" : " 10x " : Zbs ="10es13.5)') surf(index)%Zbs(1:Nfou) + if (Nbnf > 0) then + write (ounit, '(" : " 10x " : Bnim ="10i13 )') surf(index)%Bnim(1:NBnf) + write (ounit, '(" : " 10x " : Bnin ="10i13 )') surf(index)%Bnin(1:NBnf) + write (ounit, '(" : " 10x " : Bnc ="10es13.5)') surf(index)%Bnc(1:NBnf) + write (ounit, '(" : " 10x " : Bns ="10es13.5)') surf(index)%Bns(1:NBnf) + endif + endif + + surf(index)%Nteta = Nteta + surf(index)%Nzeta = Nzeta + + if (index == plasma) then + Nfp = surf(plasma)%Nfp + surf_Nfp = Nfp ! local surface Nfp + select case (IsSymmetric) + case (0) + surf_Nfp = 1 ! reset Nfp to 1; + symmetry = 0 + case (1) ! plasma and coil periodicity enabled; + symmetry = 0 + case (2) ! stellarator symmetry enforced; + symmetry = 1 + end select + + SALLOCATE(cosnfp, (1:Nfp), zero) + SALLOCATE(sinnfp, (1:Nfp), zero) + do ip = 1, Nfp + cosnfp(ip) = cos((ip - 1)*pi2/Nfp) + sinnfp(ip) = sin((ip - 1)*pi2/Nfp) + enddo + ! discretefactor = discretefactor/Nfp + surf(index)%Nzeta = Nzeta*surf_Nfp*2**symmetry ! the total number from [0, 2pi] + discretefactor = (pi2/surf(plasma)%Nteta)*(pi2/surf(plasma)%Nzeta) + endif + + SALLOCATE(surf(index)%xx, (0:Nteta - 1, 0:Nzeta - 1), zero) !x coordinates; + SALLOCATE(surf(index)%yy, (0:Nteta - 1, 0:Nzeta - 1), zero) !y coordinates + SALLOCATE(surf(index)%zz, (0:Nteta - 1, 0:Nzeta - 1), zero) !z coordinates + SALLOCATE(surf(index)%nx, (0:Nteta - 1, 0:Nzeta - 1), zero) !unit nx; + SALLOCATE(surf(index)%ny, (0:Nteta - 1, 0:Nzeta - 1), zero) !unit ny; + SALLOCATE(surf(index)%nz, (0:Nteta - 1, 0:Nzeta - 1), zero) !unit nz; + SALLOCATE(surf(index)%ds, (0:Nteta - 1, 0:Nzeta - 1), zero) !jacobian; + SALLOCATE(surf(index)%xt, (0:Nteta - 1, 0:Nzeta - 1), zero) !dx/dtheta; + SALLOCATE(surf(index)%yt, (0:Nteta - 1, 0:Nzeta - 1), zero) !dy/dtheta; + SALLOCATE(surf(index)%zt, (0:Nteta - 1, 0:Nzeta - 1), zero) !dz/dtheta; + SALLOCATE(surf(index)%pb, (0:Nteta - 1, 0:Nzeta - 1), zero) !target Bn; + SALLOCATE(surf(index)%xp, (0:Nteta - 1, 0:Nzeta - 1), zero) !dx/dzeta; + SALLOCATE(surf(index)%yp, (0:Nteta - 1, 0:Nzeta - 1), zero) !dy/dzeta; + SALLOCATE(surf(index)%zp, (0:Nteta - 1, 0:Nzeta - 1), zero) !dz/dzeta; + + surf(index)%vol = zero ! volume enclosed by plasma boundary + surf(index)%area = zero ! surface area + + ! The center point value was used to discretize grid; + do ii = 0, Nteta - 1 + teta = (ii + half)*pi2/surf(index)%Nteta + do jj = 0, Nzeta - 1 + zeta = (jj + half)*pi2/surf(index)%Nzeta + RR(0:2) = zero; ZZ(0:2) = zero + do imn = 1, surf(index)%Nfou + arg = surf(index)%bim(imn)*teta - surf(index)%bin(imn)*zeta + RR(0) = RR(0) + surf(index)%Rbc(imn)*cos(arg) + surf(index)%Rbs(imn)*sin(arg) + ZZ(0) = ZZ(0) + surf(index)%Zbc(imn)*cos(arg) + surf(index)%Zbs(imn)*sin(arg) + RR(1) = RR(1) + (-surf(index)%Rbc(imn)*sin(arg) + surf(index)%Rbs(imn)*cos(arg))*surf(index)%bim(imn) + ZZ(1) = ZZ(1) + (-surf(index)%Zbc(imn)*sin(arg) + surf(index)%Zbs(imn)*cos(arg))*surf(index)%bim(imn) + RR(2) = RR(2) - (-surf(index)%Rbc(imn)*sin(arg) + surf(index)%Rbs(imn)*cos(arg))*surf(index)%bin(imn) + ZZ(2) = ZZ(2) - (-surf(index)%Zbc(imn)*sin(arg) + surf(index)%Zbs(imn)*cos(arg))*surf(index)%bin(imn) + enddo ! end of do imn; 30 Oct 15; + szeta = sin(zeta) + czeta = cos(zeta) + xx(1:3) = (/RR(0)*czeta, RR(0)*szeta, ZZ(0)/) + xt(1:3) = (/RR(1)*czeta, RR(1)*szeta, ZZ(1)/) + xz(1:3) = (/RR(2)*czeta, RR(2)*szeta, ZZ(2)/) & + + (/-RR(0)*szeta, RR(0)*czeta, zero/) + ! minus sign for theta counterclockwise direction; + ds(1:3) = -(/xt(2)*xz(3) - xt(3)*xz(2), & + xt(3)*xz(1) - xt(1)*xz(3), & + xt(1)*xz(2) - xt(2)*xz(1)/) + dd = sqrt(sum(ds(1:3)*ds(1:3))) + ! x, y, z coordinates for the surface; + surf(index)%xx(ii, jj) = xx(1) + surf(index)%yy(ii, jj) = xx(2) + surf(index)%zz(ii, jj) = xx(3) + ! dx/dt, dy/dt, dz/dt (dt for d theta) + surf(index)%xt(ii, jj) = xt(1) + surf(index)%yt(ii, jj) = xt(2) + surf(index)%zt(ii, jj) = xt(3) + ! dx/dp, dy/dp, dz/dp (dp for d zeta(phi)) + surf(index)%xp(ii, jj) = xz(1) + surf(index)%yp(ii, jj) = xz(2) + surf(index)%zp(ii, jj) = xz(3) + ! surface normal vectors and ds for the jacobian; + surf(index)%nx(ii, jj) = ds(1)/dd + surf(index)%ny(ii, jj) = ds(2)/dd + surf(index)%nz(ii, jj) = ds(3)/dd + surf(index)%ds(ii, jj) = dd + ! using Gauss theorom; V = \int_S x \cdot n dt dz + surf(index)%vol = surf(index)%vol + surf(index)%xx(ii, jj)*ds(1) & + & + surf(index)%yy(ii, jj)*ds(2) + surf(index)%zz(ii, jj)*ds(3) + ! surface area + surf(index)%area = surf(index)%area + surf(index)%ds(ii, jj) + enddo ! end of do jj; 14 Apr 16; + enddo ! end of do ii; 14 Apr 16; + + ! print volume and area + surf(index)%vol = abs(surf(index)%vol)/3*(pi2/surf(index)%Nteta)*(pi2/surf(index)%Nzeta) + surf(index)%area = abs(surf(index)%area)*(pi2/surf(index)%Nteta)*(pi2/surf(index)%Nzeta) + if (index == plasma) then + surf(index)%vol = surf(index)%vol*surf_Nfp*2**symmetry + surf(index)%area = surf(index)%area*surf_Nfp*2**symmetry + endif + + if (myid == 0 .and. IsQuiet <= 0) then + write (ounit, '(8X": Enclosed total surface volume ="ES12.5" m^3 ; area ="ES12.5" m^2." )') & + surf(index)%vol, surf(index)%area + endif + + ! check theta direction for the plasma surface and determine the toroidal flux sign + if (index == plasma) then + dz = surf(plasma)%zz(1, 0) - surf(plasma)%zz(0, 0) + if (dz > 0) then + ! counter-clockwise + if (myid == 0) write (ounit, '(8X": The theta angle used is counter-clockwise.")') + tflux_sign = -1 + else + ! clockwise + if (myid == 0) write (ounit, '(8X": The theta angle used is clockwise.")') + tflux_sign = 1 + endif + endif + + !calculate target Bn with input harmonics; 05 Jan 17; + if (surf(index)%NBnf > 0) then + do jj = 0, Nzeta - 1 + zeta = (jj + half)*pi2/surf(index)%Nzeta + do ii = 0, Nteta - 1 + teta = (ii + half)*pi2/surf(index)%Nteta + do imn = 1, surf(index)%NBnf + arg = surf(index)%Bnim(imn)*teta - surf(index)%Bnin(imn)*zeta + surf(index)%pb(ii, jj) = surf(index)%pb(ii, jj) + surf(index)%Bnc(imn)*cos(arg) + surf(index)%Bns(imn)*sin(arg) + enddo + enddo + enddo + endif + + return + end subroutine fousurf !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! -subroutine surfcoord( index, theta, zeta, r, z) - use globals, only: dp, zero, surf - use mpi - implicit none +subroutine surfcoord(index, theta, zeta, r, z) + use globals, only: dp, zero, surf + use mpi + implicit none - INTEGER, INTENT(in) :: index - REAL, INTENT(in ) :: theta, zeta - REAL, INTENT(out) :: r, z + INTEGER, INTENT(in) :: index + REAL, INTENT(in) :: theta, zeta + REAL, INTENT(out) :: r, z - INTEGER :: imn - REAL :: arg - !-------------calculate r, z coodinates for theta, zeta------------------------------------------------ - if( .not. allocated(surf(index)%bim) ) STOP "please allocate surface data first!" + INTEGER :: imn + REAL :: arg + !-------------calculate r, z coodinates for theta, zeta------------------------------------------------ + if (.not. allocated(surf(index)%bim)) STOP "please allocate surface data first!" - r = zero; z = zero - do imn = 1, surf(index)%Nfou - arg = surf(index)%bim(imn) * theta - surf(index)%bin(imn) * zeta - R = R + surf(index)%Rbc(imn) * cos(arg) + surf(index)%Rbs(imn) * sin(arg) - Z = Z + surf(index)%Zbc(imn) * cos(arg) + surf(index)%Zbs(imn) * sin(arg) - enddo + r = zero; z = zero + do imn = 1, surf(index)%Nfou + arg = surf(index)%bim(imn)*theta - surf(index)%bin(imn)*zeta + R = R + surf(index)%Rbc(imn)*cos(arg) + surf(index)%Rbs(imn)*sin(arg) + Z = Z + surf(index)%Zbc(imn)*cos(arg) + surf(index)%Zbs(imn)*sin(arg) + enddo - return + return end subroutine surfcoord !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! diff --git a/sources/saving.f90 b/sources/saving.f90 index 85fcbca..d929bd9 100644 --- a/sources/saving.f90 +++ b/sources/saving.f90 @@ -26,7 +26,7 @@ subroutine saving !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! - INTEGER :: ii, jj, icoil, NF, ip, is, cs, Npc + INTEGER :: ii, jj, icoil, NF,NCP, ip, is, cs, Npc ! the following are used by the macros HWRITEXX below; do not alter/remove; INTEGER :: hdfier, rank @@ -55,6 +55,7 @@ subroutine saving if(allocated(t1C)) deriv(1:Ndof,7) = t1C(1:Ndof) if(allocated(t1T)) deriv(1:Ndof,8) = t1T(1:Ndof) if(allocated(t1N)) deriv(1:Ndof,9) = t1N(1:Ndof) + if(allocated(t1Str)) deriv(1:Ndof,10)=t1Str(1:Ndof) endif !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! @@ -93,6 +94,7 @@ subroutine saving HWRITEIV( 1 , case_bnormal , case_bnormal ) HWRITEIV( 1 , case_length , case_length ) HWRITEIV( 1 , case_curv , case_curv ) + HWRITEIV( 1 , case_straight , case_straight ) HWRITEIV( 1 , curv_alpha , curv_alpha ) HWRITERV( 1 , weight_bnorm , weight_bnorm ) HWRITERV( 1 , weight_bharm , weight_bharm ) @@ -107,6 +109,7 @@ subroutine saving HWRITERV( 1 , weight_inorm , weight_inorm ) HWRITERV( 1 , weight_mnorm , weight_mnorm ) HWRITERV( 1 , weight_curv , weight_curv ) + HWRITERV( 1 , weight_straight, weight_straight ) HWRITERV( 1 , weight_ccsep , weight_ccsep ) HWRITERV( 1 , weight_tors , weight_tors ) HWRITERV( 1 , ccsep_alpha , ccsep_alpha ) @@ -174,12 +177,12 @@ subroutine saving HWRITERV( 1 , Gnorm , Gnorm ) HWRITERV( 1 , Mnorm , Mnorm ) HWRITERV( 1 , overlap , overlap ) - HWRITERA( iout, 12 , evolution , evolution(1:iout, 0:11) ) + HWRITERA( iout, 13 , evolution , evolution(1:iout, 0:12) ) HWRITERA( iout, Tdof , coilspace , coilspace(1:iout, 1:Tdof) ) if (allocated(deriv)) then !HWRITERA( Ndof, 6 , deriv , deriv(1:Ndof, 0:6) ) - HWRITERA( Ndof, 9 , deriv , deriv(1:Ndof, 0:9) ) + HWRITERA( Ndof, 10 , deriv , deriv(1:Ndof, 0:10) ) endif if (allocated(Bmnc)) then @@ -210,6 +213,8 @@ subroutine saving HWRITEIV( 1 , mcssep , mcssep ) HWRITEIV( 1 , icurv , icurv ) HWRITEIV( 1 , mcurv , mcurv ) + HWRITEIV( 1 , istr , istr ) + HWRITEIV( 1 , mstr , mstr ) HWRITEIV( 1 , iccsep , iccsep ) HWRITEIV( 1 , mccsep , mccsep ) HWRITEIV( 1 , itors , itors ) @@ -238,6 +243,40 @@ subroutine saving HWRITERV( 1 , time_optimize , time_optimize ) HWRITERV( 1 , time_postproc , time_postproc ) + ! do icoil = 1,Ncoils + ! select case(icoil) + ! case(1) + ! HWRITERV( coil(icoil)%NS , curvature_1, coil(icoil)%curvature(1:coil(icoil)%NS ) ) + ! HWRITERV( coil(icoil)%NS , straight_1, coil(icoil)%straight (1:coil(icoil)%NS ) ) + ! case(2) + ! HWRITERV( coil(icoil)%NS , curvature_2, coil(icoil)%curvature(1:coil(icoil)%NS ) ) + ! HWRITERV( coil(icoil)%NS , straight_2, coil(icoil)%straight (1:coil(icoil)%NS ) ) + ! case(3) + ! HWRITERV( coil(icoil)%NS , curvature_3, coil(icoil)%curvature(1:coil(icoil)%NS ) ) + ! HWRITERV( coil(icoil)%NS , straight_3, coil(icoil)%straight (1:coil(icoil)%NS ) ) + ! case(4) + ! HWRITERV( coil(icoil)%NS , curvature_4, coil(icoil)%curvature(1:coil(icoil)%NS ) ) + ! HWRITERV( coil(icoil)%NS , straight_4, coil(icoil)%straight (1:coil(icoil)%NS ) ) + ! case(5) + ! HWRITERV( coil(icoil)%NS , curvature_5, coil(icoil)%curvature(1:coil(icoil)%NS ) ) + ! HWRITERV( coil(icoil)%NS , straight_5, coil(icoil)%straight (1:coil(icoil)%NS ) ) + ! case(6) + ! HWRITERV( coil(icoil)%NS , curvature_6, coil(icoil)%curvature(1:coil(icoil)%NS ) ) + ! HWRITERV( coil(icoil)%NS , straight_6, coil(icoil)%straight (1:coil(icoil)%NS ) ) + ! case(7) + ! HWRITERV( coil(icoil)%NS , curvature_7, coil(icoil)%curvature(1:coil(icoil)%NS ) ) + ! HWRITERV( coil(icoil)%NS , straight_7, coil(icoil)%straight (1:coil(icoil)%NS ) ) + ! case(8) + ! HWRITERV( coil(icoil)%NS , curvature_8, coil(icoil)%curvature(1:coil(icoil)%NS ) ) + ! HWRITERV( coil(icoil)%NS , straight_8, coil(icoil)%straight (1:coil(icoil)%NS ) ) + ! case(9) + ! HWRITERV( coil(icoil)%NS , curvature_9, coil(icoil)%curvature(1:coil(icoil)%NS ) ) + ! HWRITERV( coil(icoil)%NS , straight_9, coil(icoil)%straight (1:coil(icoil)%NS ) ) + ! end select + ! end do + + + call h5fclose_f( file_id, hdfier ) ! terminate access; FATAL( restart, hdfier.ne.0, error calling h5fclose_f ) @@ -282,6 +321,21 @@ subroutine saving write(wunit, *) "# Ic I Lc Bz (Ic control I; Lc control Bz)" write(wunit,'(I3, ES23.15, I3, ES23.15)') coil(icoil)%Ic, coil(icoil)%I, & coil(icoil)%Lc, coil(icoil)%Bz + case (coil_type_spline) + write(wunit, '(4(A6, A15, 8X))') " #Nseg", "current", "Ifree", "Length", "Lfree", "target_length"!, "k0" + !write(wunit,'(2X, I4, ES23.15, 3X, I3, ES23.15, 3X, I3, ES23.15, ES23.15)') & + !coil(icoil)%NS, coil(icoil)%I, coil(icoil)%Ic, coil(icoil)%L, coil(icoil)%Lc, coil(icoil)%Lo, coil(icoil)%k0 + write(wunit,'(2X, I4, ES23.15, 3X, I3, ES23.15, 3X, I3, ES23.15)') & + coil(icoil)%NS, coil(icoil)%I, coil(icoil)%Ic, coil(icoil)%L, coil(icoil)%Lc, coil(icoil)%Lo !,coil(icoil)%k0 + NCP = Splines(icoil)%NCP ! shorthand; + write(wunit, *) "#NCP " + write(wunit, '(I3,I3)') NCP + write(wunit, *) "knot vector" + write(wunit, 1000) Splines(icoil)%vect + write(wunit, *) "#Control points for coils ( x;y;z) " + write(wunit, 1000) Splines(icoil)%Cpoints(0:NCP-1) + write(wunit, 1000) Splines(icoil)%Cpoints(NCP:2*NCP-1) + write(wunit, 1000) Splines(icoil)%Cpoints(NCP*2:3*NCP-1) case default FATAL(restart, .true., not supported coil types) end select @@ -300,7 +354,7 @@ subroutine saving write(funit,'("mirror NIL")') do icoil = 1, Ncoils ! will only write x,y,z in cartesian coordinates - if (coil(icoil)%type /= 1) cycle + if ((coil(icoil)%type /= 1) .AND. (coil(icoil)%type /= coil_type_spline)) cycle ! check if the coil is stellarator symmetric select case (coil(icoil)%symm) case ( 0 ) diff --git a/sources/solvers.f90 b/sources/solvers.f90 index 9c7c37c..c7cf006 100644 --- a/sources/solvers.f90 +++ b/sources/solvers.f90 @@ -96,6 +96,12 @@ subroutine solvers write(ounit, '(8X,": ccsep_alpha = ", ES12.5, "; ccsep_beta = ", ES12.5, & & "; r_delta = ", ES12.5)') ccsep_alpha, ccsep_beta, r_delta endif + ! straight-out + if (weight_straight > machprec) then + write(ounit, '(8X,": weight_straight = ", ES12.5, "; case_straight = ", I1)') weight_straight, case_straight + write(ounit, '(8X,": str_alpha = ", ES12.5)') str_alpha + write(ounit, '(8X,": str_k0 = ", ES12.5)') str_k0 + endif endif if (abs(case_optimize) >= 1) call AllocData(1) @@ -116,46 +122,49 @@ subroutine solvers write(ounit, *) "------------- Iterations ------------------------" endif - if (lim_out .eq. 0) then - if (myid == 0) then + if (myid == 0) then + if (lim_out .eq. 0) then if (IsQuiet < 0) then - write(ounit, '("output : "A6" : "12(A12," ; "))') "iout", "mark", "chi", "dE_norm", & - "Bnormal", "Bmn harmonics", "tor. flux", "coil length", "c-s sep.", "curvature", & - "c-c sep.", "torsion", "nissin" + write(ounit, '("output : "A6" : "13(A12," ; "))') "iout", "mark", "chi", "dE_norm", & + "Bnormal", "Bn_mn harm.", "tor. flux", "coil length", "c-s sep.", "curvature", & + "c-c sep.", "torsion", "nissin", "straight" else write(ounit, '("output : "A6" : "4(A15," ; "))') "iout", "mark", "total function", "||gradient||", "Bnormal" endif - endif - else - write(ounit, '("output : "A6" : "3(A12," ; "))',advance="no") "iout", "mark", "chi", "dE_norm" - if ( weight_bnorm .ne. 0.0 ) then - write(ounit, '(A12," ; ")',advance="no") "Bnormal" - endif - if ( weight_bharm .ne. 0.0 ) then - write(ounit, '(A12," ; ")',advance="no") "Bmn harmonics" - endif - if ( weight_tflux .ne. 0.0 ) then - write(ounit, '(A12," ; ")',advance="no") "tor. flux" - endif - if ( weight_ttlen .ne. 0.0 ) then - write(ounit, '(A12," ; ")',advance="no") "coil length" - endif - if ( weight_cssep .ne. 0.0 ) then - write(ounit, '(A12," ; ")',advance="no") "c-s sep." - endif - if ( weight_curv .ne. 0.0 ) then - write(ounit, '(A12," ; ")',advance="no") "curvature" - endif - if ( weight_ccsep .ne. 0.0 ) then - write(ounit, '(A12," ; ")',advance="no") "c-c sep." - endif - if ( weight_tors .ne. 0.0 ) then - write(ounit, '(A12," ; ")',advance="no") "torsion" - endif - if ( weight_nissin .ne. 0.0 ) then - write(ounit, '(A12," ; ")',advance="no") "nissin" - endif - write(*,*) + else + write(ounit, '("output : "A6" : "3(A12," ; "))',advance="no") "iout", "mark", "chi", "dE_norm" + if ( weight_bnorm .ne. 0.0 ) then + write(ounit, '(A12," ; ")',advance="no") "Bnormal" + endif + if ( weight_bharm .ne. 0.0 ) then + write(ounit, '(A12," ; ")',advance="no") "Bn_mn harm." + endif + if ( weight_tflux .ne. 0.0 ) then + write(ounit, '(A12," ; ")',advance="no") "tor. flux" + endif + if ( weight_ttlen .ne. 0.0 ) then + write(ounit, '(A12," ; ")',advance="no") "coil length" + endif + if ( weight_cssep .ne. 0.0 ) then + write(ounit, '(A12," ; ")',advance="no") "c-s sep." + endif + if ( weight_curv .ne. 0.0 ) then + write(ounit, '(A12," ; ")',advance="no") "curvature" + endif + if ( weight_ccsep .ne. 0.0 ) then + write(ounit, '(A12," ; ")',advance="no") "c-c sep." + endif + if ( weight_tors .ne. 0.0 ) then + write(ounit, '(A12," ; ")',advance="no") "torsion" + endif + if ( weight_nissin .ne. 0.0 ) then + write(ounit, '(A12," ; ")',advance="no") "nissin" + endif + if ( weight_straight .ne. 0.0 ) then + write(ounit, '(A12," ; ")',advance="no") "straight" + endif + write(*,*) + endif endif call costfun(1) @@ -242,6 +251,7 @@ subroutine costfun(ideriv) cssep , t1S, t2S, weight_cssep, & specw , t1P, t2P, weight_specw, & ccsep , t1C, t2C, weight_ccsep, & + str , t1Str,t2Str,weight_straight, & tors , t1T, t2T, weight_tors, & nissin , t1N, t2N, weight_nissin, & curv , t1K, t2K, weight_curv, MPI_COMM_FOCUS @@ -265,6 +275,7 @@ subroutine costfun(ideriv) specw = zero ccsep = zero curv = zero + str = zero tors = zero nissin = zero @@ -290,6 +301,7 @@ subroutine costfun(ideriv) call length(0) call surfsep(0) call curvature(0) + call straight(0) call coilsep(0) call torsion(0) call nisscom(0) @@ -376,7 +388,7 @@ subroutine costfun(ideriv) endif endif - + !if(myid==0) write(ounit, '("-------i=4, chi = "ES23.15)') chi ! coil curvature if (weight_curv > machprec) then @@ -388,8 +400,21 @@ subroutine costfun(ideriv) t1E = t1E + weight_curv * t1K t2E = t2E + weight_curv * t2K endif - endif + endif + !if(myid==0) write(ounit, '("-------i=5, chi = "ES23.15)') chi + ! straight-out coil + if (weight_straight > machprec) then + call straight(ideriv) + chi = chi + weight_straight * str + if ( ideriv == 1 ) then + t1E = t1E + weight_straight * t1Str + elseif ( ideriv == 2 ) then + t1E = t1E + weight_straight * t1Str + t2E = t2E + weight_straight * t2Str + endif + endif + !if(myid==0) write(ounit, '("-------i=6, chi = "ES23.15)') chi ! coil surface separation; if (weight_cssep > machprec) then @@ -426,7 +451,7 @@ subroutine costfun(ideriv) !!$ endif !!$ !!$ endif - + !if(myid==0) write(ounit, '("-------i=7, chi = "ES23.15)') chi ! coil to coil separation if (weight_ccsep > machprec) then @@ -439,7 +464,7 @@ subroutine costfun(ideriv) t2E = t2E + weight_ccsep * t2C endif endif - + !if(myid==0) write(ounit, '("-------i=8, chi = "ES23.15)') chi ! torsion if (weight_tors > machprec) then @@ -452,7 +477,7 @@ subroutine costfun(ideriv) t2E = t2E + weight_tors * t2T endif endif - + !if(myid==0) write(ounit, '("-------i=9, chi = "ES23.15)') chi ! nissin if (weight_nissin > 0.0_dp) then @@ -511,7 +536,7 @@ subroutine normweight use globals, only: dp, zero, one, machprec, ounit, myid, xdof, bnorm, bharm, tflux, ttlen, cssep, specw, ccsep, & weight_bnorm, weight_bharm, weight_tflux, weight_ttlen, weight_cssep, weight_specw, weight_ccsep, & target_tflux, psi_avg, coil, Ncoils, case_length, Bmnc, Bmns, tBmnc, tBmns, weight_curv, curv, & - weight_tors, tors, weight_nissin, nissin, MPI_COMM_FOCUS + weight_tors, tors, weight_nissin, nissin,str,weight_straight, MPI_COMM_FOCUS use mpi implicit none @@ -603,6 +628,17 @@ subroutine normweight if( myid == 0 ) write(ounit, 1000) "weight_curv", weight_curv if( myid .eq. 0 .and. weight_curv < machprec) write(ounit, '("warning : weight_curv < machine_precision, curvature will not be used.")') + + endif + !-!-!-!-!-!-!-!-!-!-curv-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! + + if( weight_straight >= machprec ) then + + call straight(0) + if (abs(str) > machprec) weight_straight = weight_straight / str + if( myid == 0 ) write(ounit, 1000) "weight_straight", weight_straight + if( myid .eq. 0 .and. weight_straight < machprec) write(ounit, '("warning : weight_straight < machine_precision, straight-out coils& + will not be used.")') endif !-!-!-!-!-!-!-!-!-!-cssep-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! @@ -650,7 +686,7 @@ subroutine normweight !-!-!-!-!-!-!-!-!-!-nissin-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! - if( weight_nissin >= 0.0_dp ) then + if( weight_nissin > 0) then call nisscom(0) if (abs(nissin) > machprec) weight_nissin = weight_nissin / nissin @@ -675,54 +711,57 @@ subroutine output (mark) coil, coilspace, FouCoil, chi, t1E, bnorm, bharm, tflux, ttlen, cssep, specw, ccsep, & evolution, xdof, DoF, exit_tol, exit_signal, sumDE, curv, tors, nissin, MPI_COMM_FOCUS, IsQuiet, & weight_bnorm, weight_bharm, weight_tflux, weight_ttlen, weight_cssep, weight_curv, weight_ccsep, & - weight_tors, weight_nissin, lim_out + weight_tors, weight_nissin, lim_out, weight_straight, str, coil_type_spline, Splines use mpi implicit none REAL, INTENT( IN ) :: mark - INTEGER :: idof, NF, icoil + INTEGER :: idof, NF, icoil, NCP iout = iout + 1 FATAL( output , iout > Nouts+2, maximum iteration reached ) - if (lim_out .eq. 0) then - if (myid == 0) then - if (IsQuiet < 0) then - write(ounit, '("output : "I6" : "12(ES12.5," ; "))') iout, mark, chi, sumdE, bnorm, bharm, & - tflux, ttlen, cssep, curv, ccsep, tors, nissin - else - write(ounit, '("output : "I6" : "4(ES15.8," ; "))') iout, mark, chi, sumdE, bnorm - endif - endif - else - write(ounit, '("output : "I6" : "3(ES12.5," ; "))',advance="no") iout, mark, chi, sumdE - if ( weight_bnorm .ne. 0.0 ) then - write(ounit, '(ES12.5," ; ")',advance="no") bnorm - endif - if ( weight_bharm .ne. 0.0 ) then - write(ounit, '(ES12.5," ; ")',advance="no") bharm - endif - if ( weight_tflux .ne. 0.0 ) then - write(ounit, '(ES12.5," ; ")',advance="no") tflux - endif - if ( weight_ttlen .ne. 0.0 ) then - write(ounit, '(ES12.5," ; ")',advance="no") ttlen - endif - if ( weight_cssep .ne. 0.0 ) then - write(ounit, '(ES12.5," ; ")',advance="no") cssep - endif - if ( weight_curv .ne. 0.0 ) then - write(ounit, '(ES12.5," ; ")',advance="no") curv - endif - if ( weight_ccsep .ne. 0.0 ) then - write(ounit, '(ES12.5," ; ")',advance="no") ccsep - endif - if ( weight_tors .ne. 0.0 ) then - write(ounit, '(ES12.5," ; ")',advance="no") tors - endif - if ( weight_nissin .ne. 0.0 ) then - write(ounit, '(ES12.5," ; ")',advance="no") nissin - endif - write(*,*) - endif + if (myid == 0) then + if (lim_out .eq. 0) then + if (IsQuiet < 0) then + write(ounit, '("output : "I6" : "13(ES12.5," ; "))') iout, mark, chi, sumdE, bnorm, bharm, & + tflux, ttlen, cssep, curv, ccsep, tors, nissin, str + else + write(ounit, '("output : "I6" : "4(ES15.8," ; "))') iout, mark, chi, sumdE, bnorm + endif + else + write(ounit, '("output : "I6" : "3(ES12.5," ; "))',advance="no") iout, mark, chi, sumdE + if ( weight_bnorm .ne. 0.0 ) then + write(ounit, '(ES12.5," ; ")',advance="no") bnorm + endif + if ( weight_bharm .ne. 0.0 ) then + write(ounit, '(ES12.5," ; ")',advance="no") bharm + endif + if ( weight_tflux .ne. 0.0 ) then + write(ounit, '(ES12.5," ; ")',advance="no") tflux + endif + if ( weight_ttlen .ne. 0.0 ) then + write(ounit, '(ES12.5," ; ")',advance="no") ttlen + endif + if ( weight_cssep .ne. 0.0 ) then + write(ounit, '(ES12.5," ; ")',advance="no") cssep + endif + if ( weight_curv .ne. 0.0 ) then + write(ounit, '(ES12.5," ; ")',advance="no") curv + endif + if ( weight_ccsep .ne. 0.0 ) then + write(ounit, '(ES12.5," ; ")',advance="no") ccsep + endif + if ( weight_tors .ne. 0.0 ) then + write(ounit, '(ES12.5," ; ")',advance="no") tors + endif + if ( weight_nissin .ne. 0.0 ) then + write(ounit, '(ES12.5," ; ")',advance="no") nissin + endif + if ( weight_straight .ne. 0.0 ) then + write(ounit, '(ES12.5," ; ")',advance="no") str + endif + write(*,*) + endif + endif ! save evolution data; if (allocated(evolution)) then @@ -738,6 +777,7 @@ subroutine output (mark) evolution(iout,9) = ccsep evolution(iout,10)= tors evolution(iout,11)= nissin + evolution(iout,12) = str endif ! exit the optimization if no obvious changes in past 5 outputs; 07/20/2017 @@ -762,6 +802,9 @@ subroutine output (mark) coilspace(iout, idof+1:idof+NF ) = FouCoil(icoil)%zs(1:NF) ; idof = idof + NF !!$ case default !!$ FATAL(output, .true., not supported coil types) + case (coil_type_spline) + NCP = Splines(icoil)%NCP + coilspace(iout, idof+1:idof+3*NCP) = Splines(icoil)%Cpoints(0:3*NCP-1) ; idof = idof + 3*NCP end select enddo !!$ FATAL( output , idof .ne. Tdof, counting error in restart ) diff --git a/sources/splines.f90 b/sources/splines.f90 new file mode 100644 index 0000000..532a660 --- /dev/null +++ b/sources/splines.f90 @@ -0,0 +1,509 @@ +!title (boundary) ! Spline representation + +!latex \briefly{Calculates the basis functions and their derivatives for the spline representation of the coils.} + +!latex \calledby{\link{dataalloc}} +!latex \calls{\link{}} + +!latex \subsection{Overview} +!latex \item[1.] eval_spline_basis computes the basis functions of third order using the Cox-de Boor algorithm. + +!latex\[ +!latexB_{i,k}(x) = \frac{x-t_i}{t_{i+k}-t_i}B_{i,k-1}(x) + \frac{t_{i+k+1} - x}{t_{i+k+1}-t_{i+1}}B_{i+1,k-1}(x) +!latex\] +!latex\[ +!latexB_{i,0}(x) = \begin{cases} 1, \hspace{2em} t_i < x < t_{i+1} \\ 0, \hspace{2em} otherwise \end{cases} +!latex\] + +!latex It also stores the basis functions of first and second order as they are used in the computation of the derivatives +!latex \item[2.] eva_spline_basis1 computes the first derivatives of the basis functions obtained analitically as + +!latex \begin{equation*} +!latex \hspace{-4cm} \scalebox{0.6}{ $ +!latex \frac{\partial B_{i,3}(x)}{\partial x} = \begin{cases} \frac{1}{t_{i+3}-t_i}B_{i,2}(x) + \frac{x-t_i}{t_{i+3}-t_i}\left[\frac{1}{t_{i+2}-t_i}B_{i,1}(x) + \frac{x-t_i}{t_{i+2}-t_i}\frac{1}{t_{i+1}-t_i} \right], \hspace{56.5em} t_i < x < t_{i+1} \\ +!latex \frac{1}{t_{i+3}-t_i}B_{i,2}(x) + \frac{x-t_i}{t_{i+3}-t_i}\left[\frac{1}{t_{i+2}-t_i}B_{i,1}(x) - \frac{x-t_i}{t_{i+2}-t_i}\frac{1}{t_{i+2}-t_{i+1}} - \frac{1}{t_{i+3}-t_{i+1}}B_{i+1,1}(x) + \frac{t_{i+3}-x}{t_{i+3}-t_{i+1}}\frac{1}{t_{i+2}-t_{i+1}} \right] - \frac{1}{t_{i+4}-t_{i+1}}B_{i+1,2}(x) + \frac{t_{i+4}-x}{t_{i+4}-t_{i+1}}\left[\frac{1}{t_{i+3}-t_{i+1}}B_{i+1,1}(x) + \frac{x - t_{i+1}}{t_{i+3}-t_{i+1}}\frac{1}{t_{i+2}-t_{i+1}} \right], \hspace{5.5em} t_{i+1} < x < t_{i+2} \\ +!latex \frac{1}{t_{i+3}-t_i}B_{i,2}(x) + \frac{x-t_i}{t_{i+3}-t_i}\left[-\frac{1}{t_{i+3}-t_{i+1}}B_{i+1,1}(x) - \frac{t_{i+3}-x}{t_{i+3}-t_{i+1}}\frac{1}{t_{i+3}-t_{i+2}} \right] - \frac{1}{t_{i+4}-t_{i+1}}B_{i+1,2}(x) + \frac{t_{i+4}-x}{t_{i+4}-t_{i+1}}\left[\frac{1}{t_{i+3}-t_{i+1}}B_{i+1,1}(x) - \frac{x-t_{i+1}}{t_{i+3}-t_{i+1}}\frac{1}{t_{i+3}-t_{i+2}} - \frac{1}{t_{i+4}-t_{i+2}}B_{i+2,1}(x) + \frac{t_{i+4}-x}{t_{i+4}-t_{i+2}}\frac{1}{t_{i+3}-t_{i+2}} \right], \hspace{2em} t_{i+2} < x < t_{i+3} \\ +!latex -\frac{1}{t_{i+4}-t_{i+1}}B_{i+1,2}(x) + \frac{t_{i+4}-x}{t_{i+4}-t_{i+1}}\left[-\frac{1}{t_{i+4}-t_{i+2}}B_{i+2,1}(x) - \frac{t_{i+4}-x}{t_{i+4}-t_{i+2}}\frac{1}{t_{i+4}-t_{i+3}} \right], \hspace{48.5em} t_{i+3} < x < t_{i+4} \\ +!latex \end{cases}$ +!latex } +!latex \end{equation*} + +!latex \item[3.] eva_spline_basis2 computes the second derivatives of the basis functions obtained analitically as + +!latex \begin{equation*} +!latex \hspace{-4cm} \scalebox{0.5}{ $ +!latex \frac{\partial^2 B_{i,3}(x)}{\partial x^2} = \begin{cases} \frac{2}{t_{i+3}-t_i}\left[\frac{1}{t_{i+2}-t_i}B_{i,1}(x) + \frac{x-t_i}{t_{i+2}-t_i}\frac{1}{t_{i+1}-t_i} \right] +\frac{x-t_i}{t_{i+3}-t_i}\frac{2}{t_{i+1}-t_i}\frac{1}{t_{i+2}-t_i} , \hspace{71em} t_i < x < t_{i+1} \\ +!latex \frac{2}{t_{i+3}-t_i}\left[\frac{1}{t_{i+2}-t_i}B_{i,1}(x) - \frac{x-t_i}{t_{i+2}-t_i}\frac{1}{t_{i+2}-t_{i+1}} - \frac{1}{t_{i+3}-t_{i+1}}B_{i+1,1}(x) + \frac{t_{i+3}-x}{t_{i+3}-t_{i+1}}\frac{1}{t_{i+2}-t_{i+1}} \right] - \frac{x-t_i}{t_{i+3}-t_i}\frac{1}{t_{i+2}-t_{i+1}}\left[ \frac{2}{t_{i+2}-t_i} + \frac{2}{t_{i+3}-t_{i+1}}\right] - \frac{2}{t_{i+4}-t_{i+1}}\left[\frac{1}{t_{i+3}-t_{i+1}}B_{i+1,1}(x) + \frac{x-t_{i+1}}{t_{i+3}-t_{i+1}}\frac{1}{t_{i+2}-t_{i+1}} \right] + \frac{t_{i+4}-x}{t_{i+4}-t_{i+1}}\frac{2}{t_{i+3}-t_{i+1}}\frac{1}{t_{i+2}-t_{i+1}}, \hspace{9.5em} t_{i+1} < x < t_{i+2} \\ +!latex \frac{2}{t_{i+3}-t_i}\left[-\frac{1}{t_{i+3}-t_{i+1}}B_{i+1,1}(x) - \frac{t_{i+3}-x}{t_{i+3}-t_{i+1}}\frac{1}{t_{i+3}-t_{i+2}} \right] + \frac{x-t_i}{t_{i+3}-t_i}\frac{2}{t_{i+3}-t_{i+1}}\frac{1}{t_{i+3}-t_{i+2}} - \frac{2}{t_{i+4}-t_{i+1}}\left[\frac{1}{t_{i+3}-t_{i+1}}B_{i+1,1}(x) - \frac{x-t_{i+1}}{t_{i+3}-t_{i+1}}\frac{1}{t_{i+3}-t_{i+2}} - \frac{1}{t_{i+4}-t_{i+2}}B_{i+2,1}(x) + \frac{t_{i+4}-x}{t_{i+4}-t_{i+2}}\frac{1}{t_{i+3}-t_{i+2}} \right] - \frac{t_{i+4}-x}{t_{i+4}-t_{i+1}}\frac{1}{t_{i+3}-t_{i+2}}\left[\frac{2}{t_{i+3}-t_{i+1}} + \frac{2}{t_{i+4}-t_{i+2}} \right], \hspace{5em} t_{i+2} < x < t_{i+3} \\ +!latex -\frac{2}{t_{i+4}-t_{i+1}}\left[-\frac{1}{t_{i+4}-t_{i+2}}B_{i+2,1}(x) - \frac{t_{i+4}-x}{t_{i+4}-t_{i+2}}\frac{1}{t_{i+4}-t_{i+3}} \right] + \frac{t_{i+4}-x}{t_{i+4}-t_{i+1}}\frac{2}{t_{i+4}-t_{i+2}}\frac{1}{t_{i+4}-t_{i+3}}, \hspace{62.5em} t_{i+3} < x < t_{i+4} \\ +!latex \end{cases}$ +!latex } +!latex \end{equation*} + +!latex \item[4.] enforce_spline_periodicity ensures the periodicity of the spline by making sure the first three control points are used in place of the last three during optimization + +SUBROUTINE eval_spline_basis(icoil) +! Compute the basis functions of order 0,1,2 and 3 for the values of t contained in eval_points and store the value in Splines + use globals, only : dp, zero, one, myid, ounit, sqrtmachprec, IsQuiet,coil,Splines,astat + implicit none + integer :: N,i,j,icoil + REAL,allocatable :: eval_points(:) + integer :: icp,ipos + + SALLOCATE(eval_points, (0:coil(icoil)%NS-1),zero) + + N = coil(icoil)%NS + eval_points = Splines(icoil)%eval_points + + Splines(icoil)%basis_0 = 0 + Splines(icoil)%basis_1 = 0 + Splines(icoil)%basis_2 = 0 + Splines(icoil)%basis_3 = 0 + + do ipos=0,N-1 + do icp=0,Splines(icoil)%NCP+2 + if(eval_points(ipos)>=Splines(icoil)%vect(icp).AND.eval_points(ipos)=vect(icp).AND.eval_points(ipos)=vect(icp+1).AND.eval_points(ipos)=vect(icp+2).AND.eval_points(ipos)=vect(icp+3).AND.eval_points(ipos)=vect(icp).AND.eval_points(ipos)=vect(icp+1).AND.eval_points(ipos)=vect(icp+2).AND.eval_points(ipos)=vect(icp+3).AND.eval_points(ipos) 0 .AND. & + ABS((Splines(icoil)%db_dt(ipos,icp) - (bas_temp_3(ipos,icp) - Splines(icoil)%basis_3(ipos,icp))/2.0/sigma ))> 0 ) & + write (ounit,'("first pos " I7.3 " CP " I7.3 " basis 3 der" E15.7 ",rel " E15.7 " " E15.7 " " E15.7)') ipos,icp, & + ( Splines(icoil)%db_dt(ipos,icp) - (bas_temp_3(ipos,icp) - Splines(icoil)%basis_3(ipos,icp))/2.0/sigma ), & + ((Splines(icoil)%db_dt(ipos,icp) - (bas_temp_3(ipos,icp) - Splines(icoil)%basis_3(ipos,icp))/2.0/sigma ))/ & + Splines(icoil)%db_dt(ipos,icp), & + Splines(icoil)%db_dt(ipos,icp), & + (bas_temp_3(ipos,icp) - Splines(icoil)%basis_3(ipos,icp))/2.0/sigma + endif + enddo + enddo + + + Splines(icoil)%eval_points = eval_points + Splines(icoil)%eval_points= Splines(icoil)%eval_points + eval_points*sigma + call eval_spline_basis(icoil) + call eval_spline_basis1(icoil) + db_dt_temp = Splines(icoil)%db_dt + + Splines(icoil)%eval_points = Splines(icoil)%eval_points - eval_points*2.0*sigma + call eval_spline_basis(icoil) + call eval_spline_basis1(icoil) + + + do ipos = 0,coil(icoil)%NS-1 + do icp = 0,Splines(icoil)%NCP-1 + if (myid == 0) then + if (eval_points(ipos) >0 .AND. Splines(icoil)%db_dt_2(ipos,icp) > 0 .AND. & + ABS((Splines(icoil)%db_dt_2(ipos,icp) - (db_dt_temp(ipos,icp) - Splines(icoil)%db_dt(ipos,icp))/2.0/eval_points(ipos)/sigma ))> 0 ) & + write (ounit,'("second pos " I7.3 " CP " I7.3 " basis 3 der" E15.7 ",rel " E15.7 " " E15.7 " " E15.7)') ipos,icp, & + ( Splines(icoil)%db_dt_2(ipos,icp) - (db_dt_temp(ipos,icp) - Splines(icoil)%db_dt(ipos,icp))/2.0/ eval_points(ipos)/sigma ), & + ((Splines(icoil)%db_dt_2(ipos,icp) - (db_dt_temp(ipos,icp) - Splines(icoil)%db_dt(ipos,icp))/2.0/eval_points(ipos)/sigma ))/ & + Splines(icoil)%db_dt_2(ipos,icp), & + Splines(icoil)%db_dt_2(ipos,icp), & + (db_dt_temp(ipos,icp) - Splines(icoil)%db_dt(ipos,icp))/2.0/eval_points(ipos)/sigma + endif + enddo + enddo + + Splines(icoil)%eval_points = eval_points + call eval_spline_basis(icoil) + call eval_spline_basis1(icoil) +END SUBROUTINE check_eval_spline_basis + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +SUBROUTINE check_xt_xa(icoil) + !check numerically tangent vectors + use globals, only : dp, zero, one, myid, ounit, sqrtmachprec, IsQuiet,coil,Splines,astat + implicit none + integer :: N,i,j,icoil + REAL :: eval_points(0:coil(icoil)%NS-1) + integer :: icp,ipos + REAL :: vect(0:Splines(icoil)%NT-1) + REAL :: sigma = 1E-6 + REAL :: x_temp(0:coil(icoil)%NS) + REAL :: xt_temp(0:coil(icoil)%NS) + REAL :: y_temp(0:coil(icoil)%NS) + REAL :: yt_temp(0:coil(icoil)%NS) + REAL :: z_temp(0:coil(icoil)%NS) + REAL :: zt_temp(0:coil(icoil)%NS) + + vect = Splines(icoil)%vect + N = coil(icoil)%NS + + eval_points = Splines(icoil)%eval_points + + Splines(icoil)%eval_points= Splines(icoil)%eval_points + eval_points*sigma + call eval_spline_basis(icoil) + call discoil(1) + x_temp = coil(icoil)%xx + + Splines(icoil)%eval_points= Splines(icoil)%eval_points - 2.0*eval_points*sigma + call eval_spline_basis(icoil) + call discoil(1) + + do ipos = 0,coil(icoil)%NS-1 + if (.true.) then + if (eval_points(ipos) >0 .AND. coil(icoil)%xt(ipos) > 0 ) & + write (ounit,'("xt pos " I7.3 " xt der" E15.7 ",rel " E15.7 " " E15.7 " " E15.7)') ipos, & + ( coil(icoil)%xt(ipos) - (x_temp(ipos) - coil(icoil)%xx(ipos))/2.0/ eval_points(ipos)/sigma ), & + ((coil(icoil)%xt(ipos) - (x_temp(ipos) - coil(icoil)%xx(ipos))/2.0/eval_points(ipos)/sigma ))/coil(icoil)%xt(ipos), & + coil(icoil)%xt(ipos), & + (x_temp(ipos) - coil(icoil)%xx(ipos))/2.0/eval_points(ipos)/sigma + endif + enddo + + Splines(icoil)%eval_points = eval_points + call eval_spline_basis(icoil) + call discoil(1) + + Splines(icoil)%eval_points= Splines(icoil)%eval_points + eval_points*sigma + call eval_spline_basis(icoil) + call eval_spline_basis1(icoil) + call discoil(1) + xt_temp = coil(icoil)%xt + + Splines(icoil)%eval_points= Splines(icoil)%eval_points - 2.0*eval_points*sigma + call eval_spline_basis(icoil) + call eval_spline_basis1(icoil) + call discoil(1) + + do ipos = 0,coil(icoil)%NS-1 + if (myid==0) then + if (eval_points(ipos) >0 .AND. eval_points(ipos) < 2.AND. coil(icoil)%xa(ipos) > 0 ) & + write (ounit,'("xa pos " I7.3 " xa der" E15.7 ",rel " E15.7 "an " E15.7 "num " E15.7)') ipos, & + ( coil(icoil)%xa(ipos) - (xt_temp(ipos) - coil(icoil)%xt(ipos))/2.0/eval_points(ipos)/sigma ), & + ((coil(icoil)%xa(ipos) - (xt_temp(ipos) - coil(icoil)%xt(ipos))/2.0/eval_points(ipos)/sigma ))/ & + coil(icoil)%xa(ipos), & + coil(icoil)%xa(ipos), & + (xt_temp(ipos) - coil(icoil)%xt(ipos))/2.0/eval_points(ipos)/sigma + endif + + enddo + + Splines(icoil)%eval_points = eval_points + call eval_spline_basis(icoil) + call eval_spline_basis1(icoil) + call discoil(1) + + Splines(icoil)%eval_points= Splines(icoil)%eval_points + eval_points*sigma + call eval_spline_basis(icoil) + call discoil(1) + y_temp = coil(icoil)%yy + + Splines(icoil)%eval_points= Splines(icoil)%eval_points - 2.0*eval_points*sigma + call eval_spline_basis(icoil) + call discoil(1) + + do ipos = 0,coil(icoil)%NS-1 + if (.true.) then + if (eval_points(ipos) >0 .AND. eval_points(ipos) < 2.AND. coil(icoil)%yt(ipos) > 0 ) & + write (ounit,'("yt pos " I7.3 " yt der" E15.7 ",rel " E15.7 " " E15.7 " " E15.7)') ipos, & + ( coil(icoil)%yt(ipos) - (y_temp(ipos) - coil(icoil)%yy(ipos))/2.0/eval_points(ipos)/sigma ), & + ((coil(icoil)%yt(ipos) - (y_temp(ipos) - coil(icoil)%yy(ipos))/2.0/eval_points(ipos)/sigma ))/ & + coil(icoil)%yt(ipos), & + coil(icoil)%yt(ipos), & + (y_temp(ipos) - coil(icoil)%yy(ipos))/2.0/eval_points(ipos)/sigma + endif + + enddo + + Splines(icoil)%eval_points = eval_points + call eval_spline_basis(icoil) + call discoil(1) + + Splines(icoil)%eval_points= Splines(icoil)%eval_points + eval_points*sigma + call eval_spline_basis(icoil) + call eval_spline_basis1(icoil) + call discoil(1) + yt_temp = coil(icoil)%yt + + Splines(icoil)%eval_points= Splines(icoil)%eval_points - 2.0*eval_points*sigma + call eval_spline_basis(icoil) + call eval_spline_basis1(icoil) + call discoil(1) + + do ipos = 0,coil(icoil)%NS-1 + if (myid==0) then + if (eval_points(ipos) >0 .AND. eval_points(ipos) < 2.AND. coil(icoil)%ya(ipos) > 0 ) & + write (ounit,'("ya pos " I7.3 " ya der" E15.7 ",rel " E15.7 "an " E15.7 "num " E15.7)') ipos, & + ( coil(icoil)%ya(ipos) - (yt_temp(ipos) - coil(icoil)%yt(ipos))/2.0/eval_points(ipos)/sigma ), & + ((coil(icoil)%ya(ipos) - (yt_temp(ipos) - coil(icoil)%yt(ipos))/2.0/eval_points(ipos)/sigma ))/ & + coil(icoil)%ya(ipos), & + coil(icoil)%ya(ipos), & + (yt_temp(ipos) - coil(icoil)%yt(ipos))/2.0/eval_points(ipos)/sigma + endif + + enddo + + Splines(icoil)%eval_points = eval_points + call eval_spline_basis(icoil) + call eval_spline_basis1(icoil) + call discoil(1) + + Splines(icoil)%eval_points= Splines(icoil)%eval_points + eval_points*sigma + call eval_spline_basis(icoil) + call eval_spline_basis1(icoil) + call discoil(1) + + zt_temp = coil(icoil)%zt + Splines(icoil)%eval_points= Splines(icoil)%eval_points - 2.0*eval_points*sigma + call eval_spline_basis(icoil) + call eval_spline_basis1(icoil) + call discoil(1) + + do ipos = 0,coil(icoil)%NS-1 + if (myid==0) then + if (eval_points(ipos) >0 .AND. eval_points(ipos) < 2.AND. coil(icoil)%za(ipos) > 0 ) & + write (ounit,'("za pos " I7.3 " za der" E15.7 ",rel " E15.7 "an " E15.7 "num " E15.7)') ipos, & + ( coil(icoil)%za(ipos) - (zt_temp(ipos) - coil(icoil)%zt(ipos))/2.0/eval_points(ipos)/sigma ), & + ((coil(icoil)%za(ipos) - (zt_temp(ipos) - coil(icoil)%zt(ipos))/2.0/eval_points(ipos)/sigma ))/ & + coil(icoil)%za(ipos), & + coil(icoil)%za(ipos), & + (zt_temp(ipos) - coil(icoil)%zt(ipos))/2.0/eval_points(ipos)/sigma + endif + + enddo + + Splines(icoil)%eval_points = eval_points + call eval_spline_basis(icoil) + call eval_spline_basis1(icoil) + call discoil(1) + + Splines(icoil)%eval_points= Splines(icoil)%eval_points + eval_points*sigma + call eval_spline_basis(icoil) + call discoil(1) + z_temp = coil(icoil)%zz + + Splines(icoil)%eval_points= Splines(icoil)%eval_points - 2.0*eval_points*sigma + call eval_spline_basis(icoil) + call discoil(1) + do ipos = 0,coil(icoil)%NS-1 + if (.true.) then + if (eval_points(ipos) >0 .AND. eval_points(ipos) < 2 .AND. coil(icoil)%zt(ipos) > 0 ) & + write (ounit,'("zt pos " I7.3 " zt der" E15.7 ",rel " E15.7 " " E15.7 " " E15.7)') ipos, & + ( coil(icoil)%zt(ipos) - (z_temp(ipos) - coil(icoil)%zz(ipos))/2.0/eval_points(ipos)/sigma ), & + ((coil(icoil)%zt(ipos) - (z_temp(ipos) - coil(icoil)%zz(ipos))/2.0/eval_points(ipos)/sigma ))/ & + coil(icoil)%zt(ipos), & + coil(icoil)%zt(ipos), & + (z_temp(ipos) - coil(icoil)%zz(ipos))/2.0/eval_points(ipos)/sigma + endif + + enddo + + Splines(icoil)%eval_points = eval_points + + call eval_spline_basis(icoil) + call eval_spline_basis1(icoil) + call discoil(1) + +END SUBROUTINE check_xt_xa diff --git a/sources/straight_out.f90 b/sources/straight_out.f90 new file mode 100644 index 0000000..39fde38 --- /dev/null +++ b/sources/straight_out.f90 @@ -0,0 +1,411 @@ + +!title (Straight out) ! Calculate objective functon for straight-out coils and its derivatives. (tkruger) + +!latex \briefly{The cost function will allow for straight-out coils.} + +!latex \calledby{\link{solvers}} + +!latex \section{General} +!latex Computes the straight out penalty in one of three possible forms. It is similar to the curvature penality but only applied to the outer part of the coil. +!latex \item[1.] A linear objective function +!latex \begin{eqnarray} +!latex f_{curv} = \frac{1}{N_c}\sum_{i=1}^{N_c} \int_0^{2\pi}W(t)\kappa_i dt +!latex \end{eqnarray} +!latex \item A quadratic objective function +!latex \begin{eqnarray} +!latex f_{curv} = \frac{1}{N_c}\sum_{i=1}^{N_c} \int_0^{2\pi}W(t)\kappa_i^2 dt +!latex \end{eqnarray} +!latex \item A maximum curvature objective function +!latex \begin{eqnarray} +!latex f_{curv} = \frac{1}{N_c}\sum_{i=1}^{N_c} \int_0^{2\pi} W(t) H_{\kappa_o}(\kappa_i) (\kappa_i - \kappa_o)^\alpha dt +!latex \end{eqnarray} + +!latex where $H_{\kappa_o}(\kappa_i)$ is the step function, $\kappa_o$ and $\alpha$ are user-defined parameters and $\kappa$ is the curvature at a point defined as +!latex \begin{equation} +!latex \kappa = \frac{( (z^\prime^\prime y^\prime - y^\prime^\prime z^\prime)^2 + (x^\prime^\prime z^\prime - z^\prime^\prime x^\prime)^2 + (y^\prime^\prime x^\prime - x^\prime^\prime y^\prime)^2)^\frac{1}{2}}{(x^\prime^2 + y^\prime^2 + z^\prime^2 )^\frac{3}{2}} +!latex \end{equation} +!latex In particular $\kappa_o$ has the meaning of a maximum allowed curvature and the cost function has no effect on points with a curvature smaller than this parameter. + +!latex W(t) is a weight function to enure the penalty is only applied to the outer part of the coil +!latex \begin{equation} +!latex W(t) = \begin{cases} 1, \hspace{2em} P_{xy}(t) > P_m + \beta P_d \\ 0, \hspace{2em} otherwise \end{cases} +!latex where $P_m$ is the mean squared distance of the coil projection in the x-y plane from a user-defined point and $P_d$ +!latex is a measure for half the maximum projection in the x-y plane of the coil radius. +!latex Considering the distance of the projection of a point along the coil in the x-y plane from the user-defined point ($x_0,y_0$) +!latex \begin{equation} +!latex P_{xy}(t) = (x_i(t) - x_0)^2 + (y_i(t) - y_0)^2 \ , +!latex \end{equation} +!latex they are defined as +!latex \begin{align} +!latex & P_m = \overline{P_{xy}} \ , & P_d = \frac{\max(P_{xy}) - P_m}{2} \ . +!latex \end{align} +!latex $\beta$ is a user-defined parameter that can be used to tweak the section of the coil affected by the optimization. +!latex With a value $\beta = 0$ the section of the coil with a projected distance greater than the mean value for that coil will all be affected +!latex while for a value $\beta = 2$ the new cost function will not be applied to any point. + +!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! + +! str is total penalty +! chi = chi + weight_straight*str +! t1Str is total derivative of penalty +! LM implemented +! not parallelized, at some point check to see how long takes to run +subroutine straight(ideriv) + use globals, only: dp, zero, half, pi2, machprec, ncpu, myid, ounit, MPI_COMM_FOCUS, & + coil, DoF, Ncoils, Nfixgeo, Ndof, str, t1Str, t2Str, weight_straight, FouCoil, & + mstr, istr, LM_fvec, LM_fjac,coil_type_spline,Splines,origin_surface_x,origin_surface_y,origin_surface_z,coeff_disp_straight + + implicit none + include "mpif.h" + INTEGER, INTENT(in) :: ideriv + + INTEGER :: astat, ierr, icoil, idof, ND, NF,NCP, ivec + REAL :: strAdd + + !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! + str = zero + strAdd = zero + ivec = 1 + + + + if( ideriv >= 0 ) then + + do icoil = 1, Ncoils + if( (coil(icoil)%type .ne. 1) .AND. (coil(icoil)%type .ne. coil_type_spline) ) exit ! only for Fourier + if( coil(icoil)%Lc /= 0 ) then ! if geometry is free + + call StrDeriv0(icoil,strAdd) + str = str + strAdd + + if (mstr > 0) then ! L-M format of targets + LM_fvec(istr+ivec) = weight_straight*strAdd + ivec = ivec + 1 + endif + endif + enddo + + str = str / (Ncoils - Nfixgeo + machprec) + + endif + + !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! + + if ( ideriv >= 1 ) then + + t1Str = zero + ivec = 1 + idof = 0 + icoil = 1 + !call unpacking(xdof) + !if (myid == 0 .AND. icoil==1) write(ounit,'("R**2 " 7F20.10)')(coil(icoil)%xx-origin_surface_x)**2 + & + ! (coil(icoil)%yy-origin_surface_y)**2 + do icoil = 1, Ncoils + + if( (coil(icoil)%type .ne. 1) .AND. (coil(icoil)%type .ne. coil_type_spline) ) exit ! only for Fourier + + ND = DoF(icoil)%ND + if (coil(icoil)%type == 1) NF = FouCoil(icoil)%NF + if (coil(icoil)%type == coil_type_spline) NCP = Splines(icoil)%NCP + + if ( coil(icoil)%Ic /= 0 ) then !if current is free; + idof = idof + 1 + endif + + if ( coil(icoil)%Lc /= 0 ) then !if geometry is free; + if (coil(icoil)%type == coil_type_spline) call StrDeriv1( icoil, t1Str(idof+1:idof+ND), ND, NCP ) + if (coil(icoil)%type == 1) call StrDeriv1( icoil, t1Str(idof+1:idof+ND), ND, NF ) + if (mstr > 0) then ! L-M format of targets + LM_fjac(istr+ivec, idof+1:idof+ND) = weight_straight * t1Str(idof+1:idof+ND) + ivec = ivec + 1 + endif + idof = idof + ND + endif + + enddo + FATAL( straight-out , idof .ne. Ndof, counting error in packing ) + + t1Str = t1Str / (Ncoils - Nfixgeo + machprec) + + endif + + return +end subroutine straight + +!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! +! calculate coil curvature + +subroutine StrDeriv0(icoil,strRet) + + use globals, only: dp, zero, pi2, ncpu, astat, ierr, myid, ounit, coil, NFcoil, Nseg, Ncoils, & + case_straight, str_alpha, str_k0, str_k1, str_beta, str_gamma, str_sigma, penfun_str, & + str_k1len, MPI_COMM_FOCUS,coil_type_spline,Splines, origin_surface_x, origin_surface_y, origin_surface_z,coeff_disp_straight + + + implicit none + include "mpif.h" + + INTEGER, intent(in) :: icoil + REAL , intent(out) :: strRet +!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! + INTEGER :: kseg, NS + REAL,allocatable :: strv(:), xt(:), yt(:), zt(:), xa(:), ya(:), za(:) + REAL :: hypc, str_hold, k1_use,mean_xy_distance, dispersion + + NS = coil(icoil)%NS + + SALLOCATE(strv, (0:NS-1),zero) + SALLOCATE(xt, (0:NS), zero) + SALLOCATE(yt, (0:NS), zero) + SALLOCATE(zt, (0:NS), zero) + SALLOCATE(xa, (0:NS), zero) + SALLOCATE(ya, (0:NS), zero) + SALLOCATE(za, (0:NS), zero) + xt(0:coil(icoil)%NS) = coil(icoil)%xt(0:coil(icoil)%NS) + yt(0:coil(icoil)%NS) = coil(icoil)%yt(0:coil(icoil)%NS) + zt(0:coil(icoil)%NS) = coil(icoil)%zt(0:coil(icoil)%NS) + xa(0:coil(icoil)%NS) = coil(icoil)%xa(0:coil(icoil)%NS) + ya(0:coil(icoil)%NS) = coil(icoil)%ya(0:coil(icoil)%NS) + za(0:coil(icoil)%NS) = coil(icoil)%za(0:coil(icoil)%NS) + + if ( case_straight .eq. 1 ) then + str_alpha = 0.0 + str_sigma = 1.0 + str_gamma = 1.0 + str_k1 = 0.0 + elseif ( case_straight .eq. 2 ) then + str_alpha = 0.0 + str_sigma = 1.0 + str_gamma = 2.0 + str_k1 = 0.0 + elseif (case_straight .eq. 3 ) then + str_sigma = 0.0 + else + ! Do nothing + endif + + FATAL( StrDeriv0, icoil .lt. 1 .or. icoil .gt. Ncoils, icoil not in right range ) + FATAL( StrDeriv0, penfun_str .ne. 1 .and. penfun_str .ne. 2 , invalid choice of penfun_str, pick 1 or 2 ) + FATAL( StrDeriv0, str_k0 < 0.0 , str_k0 cannot be negative ) + FATAL( StrDeriv0, str_k1 < 0.0 , str_k1 cannot be negative ) + FATAL( StrDeriv0, str_alpha < 0.0 , str_alpha cannot be negative ) + FATAL( StrDeriv0, str_beta < 2.0 , str_beta needs to be >= 2 ) + FATAL( StrDeriv0, str_gamma < 1.0 , str_gamma needs to be >= 1 ) + FATAL( StrDeriv0, str_sigma < 0.0 , str_sigma cannot be negative ) + FATAL( StrDeriv0, str_gamma .eq. 1.0 .and. str_k1 .ne. 0.0 , if str_gamma = 1, str_k1 must = 0 ) + + strv = zero + strRet = zero + mean_xy_distance = SUM((coil(icoil)%xx-origin_surface_x)**2 + (coil(icoil)%yy-origin_surface_y)**2)/(coil(icoil)%NS) + !dispersion = (MAXVAL((coil(icoil)%xx-origin_surface_x)**2 + (coil(icoil)%yy-origin_surface_y)**2)- & +! MINVAL((coil(icoil)%xx-origin_surface_x)**2 + (coil(icoil)%yy-origin_surface_y)**2))/2 + dispersion = (MAXVAL((coil(icoil)%xx-origin_surface_x)**2 + (coil(icoil)%yy-origin_surface_y)**2) - mean_xy_distance)/2 + do kseg = 0,coil(icoil)%NS-1 + if ((coil(icoil)%xx(kseg)-origin_surface_x)**2 + (coil(icoil)%yy(kseg)-origin_surface_y)**2 > & + (mean_xy_distance + coeff_disp_straight*dispersion)) then + strv(kseg) = sqrt( (za(kseg)*yt(kseg)-zt(kseg)*ya(kseg))**2 + (xa(kseg)*zt(kseg)-xt(kseg)*za(kseg))**2 + (ya(kseg)*xt(kseg)-yt(kseg)*xa(kseg))**2 )& + / ((xt(kseg))**2+(yt(kseg))**2+(zt(kseg))**2)**(1.5) + else + strv(kseg) = 0 + endif + enddo + !if (myid==0 .AND. icoil==1) write(ounit,'("strv "7F20.10)')strv + coil(icoil)%straight = strv + + + if ( str_k1len .eq. 1 ) then + k1_use = pi2/coil(icoil)%Lo + else + k1_use = str_k1 + endif + + do kseg = 0,NS-1 + str_hold = 0.0 + if ( str_alpha .ne. 0.0 ) then + if ( strv(kseg) > str_k0 ) then + if ( penfun_str .eq. 1 ) then + hypc = 0.5*exp( str_alpha*( strv(kseg) - str_k0 ) ) + 0.5*exp( -1.0*str_alpha*( strv(kseg) - str_k0 ) ) + str_hold = (hypc - 1.0)**2 + else + str_hold = ( str_alpha*(strv(kseg)-str_k0) )**str_beta + endif + endif + endif + if ( str_sigma .ne. 0.0 ) then + if ( strv(kseg) > k1_use ) then + str_hold = str_hold + str_sigma*( ( strv(kseg) - k1_use )**str_gamma ) + endif + + endif + strRet = strRet + str_hold*sqrt(xt(kseg)**2+yt(kseg)**2+zt(kseg)**2) + enddo + + call lenDeriv0( icoil, coil(icoil)%L ) + + if (coil(icoil)%type==1) strRet = pi2*strRet/(NS*coil(icoil)%L) + if (coil(icoil)%type==coil_type_spline) strRet = 1.0*strRet/(NS*coil(icoil)%L) + + DALLOCATE(strv) + DALLOCATE(xt) + DALLOCATE(yt) + DALLOCATE(zt) + DALLOCATE(xa) + DALLOCATE(ya) + DALLOCATE(za) + + return + +end subroutine StrDeriv0 + +!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! +subroutine StrDeriv1(icoil, derivs, ND, NC ) !Calculate all derivatives for a coil + + use globals, only: dp, zero, pi2, coil, DoF, myid, ounit, Ncoils, & + case_straight, str_alpha, str_k0, str_k1, str_beta, str_gamma & + ,str_sigma, penfun_str, str_k1len, FouCoil, MPI_COMM_FOCUS,Splines,coil_type_spline, & + origin_surface_x, origin_surface_y, origin_surface_z,xdof,coeff_disp_straight + implicit none + !include "mpif.h" + + INTEGER, intent(in) :: icoil, ND , NC !NC is actually NCP for spline coils + REAL , intent(out) :: derivs(1:1, 1:ND) +!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! + INTEGER :: kseg, astat, ierr, doff, nff, NS,n + REAL :: dl3, xt, yt, zt, xa, ya, za, f1, f2, ncc, nss, ncp, nsp, strHold, penStr, strH, leng, hypc, hyps, str_deriv, & + k1_use, rtxrax, rtxray, rtxraz,mean_xy_distance,dispersion + REAL, dimension(1:1, 0:coil(icoil)%NS-1) :: dLx, dLy, dLz + REAL :: d1L(1:1, 1:ND) + REAL,allocatable :: dxtdDoF(:,:), dytdDoF(:,:), dztdDoF(:,:), dxadDoF(:,:), dyadDoF(:,:), dzadDoF(:,:) + + FATAL( StrDeriv1, icoil .lt. 1 .or. icoil .gt. Ncoils, icoil not in right range ) + + mean_xy_distance = SUM((coil(icoil)%xx-origin_surface_x)**2 + (coil(icoil)%yy-origin_surface_y)**2)/(coil(icoil)%NS) + !dispersion = (MAXVAL((coil(icoil)%xx-origin_surface_x)**2 + (coil(icoil)%yy-origin_surface_y)**2)- & +! MINVAL((coil(icoil)%xx-origin_surface_x)**2 + (coil(icoil)%yy-origin_surface_y)**2))/2 + + dispersion = (MAXVAL((coil(icoil)%xx-origin_surface_x)**2 + (coil(icoil)%yy-origin_surface_y)**2) - mean_xy_distance)/2 + + !write(ounit, '( " dist " F20.10 " disp " F20.10)') & +! mean_xy_distance,dispersion + derivs = zero + derivs(1,1:ND) = 0.0 + NS = coil(icoil)%NS + d1L = zero + + SALLOCATE(dxtdDoF, (0:NS,1:ND), zero) + SALLOCATE(dytdDoF, (0:NS,1:ND), zero) + SALLOCATE(dztdDoF, (0:NS,1:ND), zero) + SALLOCATE(dxadDoF, (0:NS,1:ND), zero) + SALLOCATE(dyadDoF, (0:NS,1:ND), zero) + SALLOCATE(dzadDoF, (0:NS,1:ND), zero) + + select case (coil(icoil)%type) + case(1) + do n = 1,NC + dxtdDof(0:NS,n+1) = -1*FouCoil(icoil)%smt(0:NS,n) * n + dxtdDof(0:NS,n+NC+1) = FouCoil(icoil)%cmt(0:NS,n) * n + dytdDof(0:NS,n+2*NC+2) = -1*FouCoil(icoil)%smt(0:NS,n) * n + dytdDof(0:NS,n+3*NC+2) = FouCoil(icoil)%cmt(0:NS,n) * n + dztdDof(0:NS,n+4*NC+3) = -1*FouCoil(icoil)%smt(0:NS,n) * n + dztdDof(0:NS,n+5*NC+3) = FouCoil(icoil)%cmt(0:NS,n) * n + + dxadDof(0:NS,n+1) = -1*FouCoil(icoil)%cmt(0:NS,n) * n*n + dxadDof(0:NS,n+NC+1) = -1*FouCoil(icoil)%smt(0:NS,n) * n*n + dyadDof(0:NS,n+2*NC+2) = -1*FouCoil(icoil)%cmt(0:NS,n) * n*n + dyadDof(0:NS,n+3*NC+2) = -1*FouCoil(icoil)%smt(0:NS,n) * n*n + dzadDof(0:NS,n+4*NC+3) = -1*FouCoil(icoil)%cmt(0:NS,n) * n*n + dzadDof(0:NS,n+5*NC+3) = -1*FouCoil(icoil)%smt(0:NS,n) * n*n + enddo + case(coil_type_spline) + dxtdDof(0:NS,0:ND) = Splines(icoil)%db_dt(0:NS,0:ND) + dytdDof(0:NS,0:ND) = Splines(icoil)%db_dt(0:NS,0:ND) + dztdDof(0:NS,0:ND) = Splines(icoil)%db_dt(0:NS,0:ND) + + dxadDof(0:NS,0:ND) = Splines(icoil)%db_dt_2(0:NS,0:ND) + dyadDof(0:NS,0:ND) = Splines(icoil)%db_dt_2(0:NS,0:ND) + dzadDof(0:NS,0:ND) = Splines(icoil)%db_dt_2(0:NS,0:ND) + case default + FATAL( StrDeriv1, .true. , invalid coil_type option ) + end select + + + call lenDeriv0( icoil, coil(icoil)%L ) + leng = coil(icoil)%L + call lenDeriv1( icoil, d1L(1:1,1:ND), ND ) + + if ( case_straight .eq. 1 ) then + str_alpha = 0.0 + str_sigma = 1.0 + str_gamma = 1.0 + str_k1 = 0.0 + elseif ( case_straight .eq. 2 ) then + str_alpha = 0.0 + str_sigma = 1.0 + str_gamma = 2.0 + str_k1 = 0.0 + elseif (case_straight .eq. 3 ) then + str_sigma = 0.0 + else + ! Do nothing + endif + + do kseg = 0,NS-1 + + if (((coil(icoil)%xx(kseg)-origin_surface_x)**2 + (coil(icoil)%yy(kseg)-origin_surface_y)**2) > & + (mean_xy_distance + coeff_disp_straight*dispersion)) then + + xt = coil(icoil)%xt(kseg) ; yt = coil(icoil)%yt(kseg) ; zt = coil(icoil)%zt(kseg) ; + xa = coil(icoil)%xa(kseg) ; ya = coil(icoil)%ya(kseg) ; za = coil(icoil)%za(kseg) ; + rtxrax = yt*za - zt*ya + rtxray = zt*xa - xt*za + rtxraz = xt*ya - yt*xa + f1 = sqrt( (xt*ya-xa*yt)**2 + (xt*za-xa*zt)**2 + (yt*za-ya*zt)**2 ); + f2 = (xt**2+yt**2+zt**2)**(1.5); + strHold = f1/f2 + penStr = 0.0 + str_deriv = 0.0 + + if ( str_k1len .eq. 1 ) then + k1_use = pi2/coil(icoil)%Lo + else + k1_use = str_k1 + endif + if ( penfun_str .eq. 1 ) then + if ( strHold > str_k0 ) then + hypc = 0.5*exp( str_alpha*(strHold-str_k0) ) + 0.5*exp( -1.0*str_alpha*(strHold-str_k0) ) + hyps = 0.5*exp( str_alpha*(strHold-str_k0) ) - 0.5*exp( -1.0*str_alpha*(strHold-str_k0) ) + penStr = ( hypc - 1.0 )**2 + str_deriv = 2.0*str_alpha*( hypc - 1.0 )*hyps + endif + else + if ( strHold > str_k0 ) then + penStr = (str_alpha*(strHold-str_k0))**str_beta + str_deriv = str_beta*str_alpha*( (str_alpha*(strHold-str_k0))**(str_beta-1.0) ) + endif + endif + if ( strHold > k1_use ) then + str_deriv = str_deriv + str_sigma*str_gamma*( (strHold-k1_use)**(str_gamma-1.0) ) + penStr = penStr + str_sigma*( (strHold-k1_use)**str_gamma ) + endif + + derivs(1,1:ND) = derivs(1,1:ND) + sqrt(xt**2+yt**2+zt**2)*penStr*-1.0*d1L(1,1:ND)/leng + derivs(1,1:ND) = derivs(1,1:ND) + (xt*dxtdDof(kseg,1:ND)+yt*dytdDof(kseg,1:ND)+zt*dztdDof(kseg,1:ND))*(xt**2+yt**2+zt**2)**(-.5)*penStr + derivs(1,1:ND) = derivs(1,1:ND) + sqrt(xt**2+yt**2+zt**2) & + *(-3.0*(f1/f2**2)*sqrt(xt**2+yt**2+zt**2)*(xt*dxtdDof(kseg,1:ND)+yt*dytdDof(kseg,1:ND)+zt*dztdDof(kseg,1:ND)) & + + ( rtxrax*(dytdDof(kseg,1:ND)*za - dztdDof(kseg,1:ND)*ya + yt*dzadDof(kseg,1:ND) - zt*dyadDof(kseg,1:ND)) & + + rtxray*(dztdDof(kseg,1:ND)*xa - dxtdDof(kseg,1:ND)*za + zt*dxadDof(kseg,1:ND) - xt*dzadDof(kseg,1:ND)) & + + rtxraz*(dxtdDof(kseg,1:ND)*ya - dytdDof(kseg,1:ND)*xa + xt*dyadDof(kseg,1:ND) - yt*dxadDof(kseg,1:ND)) )/(f1*f2) )*str_deriv + endif + enddo + + if (coil(icoil)%type == 1) derivs = pi2*derivs/(NS*leng) + if (coil(icoil)%type == coil_type_spline )derivs = derivs/(NS*leng) + + + !if (myid==0 .AND. icoil==1) write(ounit,'(" str derivs " 7F20.10)')derivs + !if (myid==0 .AND. icoil==3) write(ounit,'(" 3 derivs " 7F20.10)')derivs + + return + +end subroutine StrDeriv1 + diff --git a/sources/surface.f90 b/sources/surface.f90 index e955351..3f9cd90 100644 --- a/sources/surface.f90 +++ b/sources/surface.f90 @@ -1,7 +1,7 @@ ! This is the overall function to handle surfaces SUBROUTINE surface use globals, only : dp, myid, ounit, machprec, surf, plasma, limiter, input_surf, limiter_surf, & - psurf, weight_cssep, MPI_COMM_FOCUS + psurf, weight_cssep, MPI_COMM_FOCUS,plasma_surf_boozer,case_surface use mpi implicit none @@ -24,10 +24,17 @@ SUBROUTINE surface ! read the plasma surface inquire( file=trim(input_surf), exist=exist) FATAL( surface, .not.exist, input_surf does not exist ) - call fousurf( input_surf, plasma ) + + if (case_surface == plasma_surf_boozer) call rdbooz( input_surf, plasma ) + if (case_surface .ne. plasma_surf_boozer) call fousurf( input_surf, plasma ) ! read the limiter surface - if (limiter /= plasma) then + if (limiter /= plasma .and. case_surface == plasma_surf_boozer) then + inquire( file=trim(limiter_surf), exist=exist) + FATAL( surface, .not.exist, limiter_surf does not exist ) + FATAL( surface, limiter <= plasma, something goes wrong the surface indexing ) + call rdbooz( limiter_surf, limiter ) + elseif (limiter /= plasma .and. case_surface .ne. plasma_surf_boozer) then inquire( file=trim(limiter_surf), exist=exist) FATAL( surface, .not.exist, limiter_surf does not exist ) FATAL( surface, limiter <= plasma, something goes wrong the surface indexing ) diff --git a/sources/surfsep.f90 b/sources/surfsep.f90 index 9f1a513..35e772a 100644 --- a/sources/surfsep.f90 +++ b/sources/surfsep.f90 @@ -73,7 +73,8 @@ SUBROUTINE surfsep(ideriv) !------------------------------------------------------------------------------------------------------ use globals, only: dp, zero, half, pi2, machprec, ncpu, myid, ounit, & coil, DoF, Ncoils, Nfixgeo, Ndof, cssep, t1S, t2S, psurf, surf, & - icssep, mcssep, LM_fvec, LM_fjac, weight_cssep, Nteta, Nzeta, Nfp, MPI_COMM_FOCUS + icssep, mcssep, LM_fvec, LM_fjac, weight_cssep, Nteta, Nzeta, Nfp, MPI_COMM_FOCUS,coil_type_spline + use mpi implicit none @@ -90,11 +91,11 @@ SUBROUTINE surfsep(ideriv) discretefactor = (pi2/Nteta) * (pi2/Nzeta) lcssep = zero totalcoil = zero - + if( ideriv >= 0 ) then ivec = 1 do icoil = 1, Ncoils - if (coil(icoil)%type /= 1) cycle ! skip for other coils + if ((coil(icoil)%type /= 1) .AND. (coil(icoil)%type /= coil_type_spline)) cycle ! skip for other coils coilsum = zero if ( coil(icoil)%Lc /= 0 ) then if ( coil(icoil)%symm == 0 ) then @@ -147,7 +148,7 @@ SUBROUTINE surfsep(ideriv) endif if ( coil(icoil)%Lc /= 0 ) then ! if geometry is free; - if (coil(icoil)%type == 1) then ! skip for other coils + if ((coil(icoil)%type == 1) .OR. (coil(icoil)%type == coil_type_spline)) then ! skip for other coils do jzeta = 0, Nzeta - 1 do iteta = 0, Nteta - 1 if( myid.ne.modulo(jzeta*Nteta+iteta,ncpu) ) cycle ! parallelization loop; @@ -176,7 +177,7 @@ SUBROUTINE surfsep(ideriv) do idof = 1, Ndof t1S(idof) = sum(jac(1:Ncoils, idof)) * discretefactor / (totalcoil + machprec) enddo - + endif return diff --git a/sources/torflux.f90 b/sources/torflux.f90 index fcad91d..6129a52 100644 --- a/sources/torflux.f90 +++ b/sources/torflux.f90 @@ -236,7 +236,7 @@ subroutine bpotential0(icoil, iteta, jzeta, tAx, tAy, tAz) ! Discretizing factor is includeed; coil(icoil)%dd(kseg) !------------------------------------------------------------------------------------------------------ use globals, only: dp, coil, surf, Ncoils, Nteta, Nzeta, MPI_COMM_FOCUS, & - zero, myid, ounit, plasma, Nfp, cosnfp, sinnfp, two, bsconstant + zero, myid, ounit, plasma, Nfp, cosnfp, sinnfp, two, bsconstant,coil_type_spline use mpi implicit none @@ -281,7 +281,7 @@ subroutine bpotential0(icoil, iteta, jzeta, tAx, tAy, tAz) zz = surf(isurf)%zz(iteta,jzeta) * (-1)**is Ax = zero; Ay = zero; Az = zero select case (coil(icoil)%type) - case(1) + case(1,coil_type_spline) ! Fourier coils do kseg = 0, coil(icoil)%NS-1 dlx = xx - coil(icoil)%xx(kseg) @@ -324,7 +324,7 @@ subroutine bpotential1(icoil, iteta, jzeta, tAx, tAy, tAz, ND) ! Discretizing factor is includeed; coil(icoil)%dd(kseg) !------------------------------------------------------------------------------------------------------ use globals, only: dp, coil, DoF, surf, NFcoil, Ncoils, Nteta, Nzeta, & - zero, myid, ounit, plasma, Nfp, cosnfp, sinnfp, MPI_COMM_FOCUS + zero, myid, ounit, plasma, Nfp, cosnfp, sinnfp, MPI_COMM_FOCUS,coil_type_spline use mpi implicit none @@ -373,7 +373,7 @@ subroutine bpotential1(icoil, iteta, jzeta, tAx, tAy, tAz, ND) zz = surf(isurf)%zz(iteta,jzeta) * (-1)**is Ax = zero; Ay = zero; Az = zero select case (coil(icoil)%type) - case(1) + case(1,coil_type_spline) ! Fourier coils do kseg = 0, NS-1 dlx = xx - coil(icoil)%xx(kseg)